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Introduction 


This text was written to accompany a series of lectures on computational methods for 
fluid dynamics problems, which I gave in the Space Science Division at NASA Ames 
Research Center in July and August of 1986. Each chapter in the text corresponds to a 
particular lecture. Sufficient interest was expressed in the subject that Lynda Haines, the 
User Services Manager of the Numerical Aerodynamic Simulator (NAS) Project, asked 
me to prepare a series of videotaped lectures to be distributed along with this text. 
Anyone interested in obtaining these tapes should contact NAS User Services at NASA 
Ames Research Center, Moffett Field, CA 94035. 

This lecture series covers the basic principles of computational fluid dynamics (ab- 
breviated as “CFD” throughout). The lectures are designed to teach an inexperienced 
person everything he or she needs to know to create a time dependent numerical model 
of fluid flow, in one or more dimensions. I say “a model” because there is no unique 
way to construct numerical models, and the number of such models which have been 
constructed greatly exceeds the number of practitioners in the field. 

As there are so many different approaches to simulating fluid flow, it is appropriate 
to consider what will and will not be covered by these lectures. Lagrangian methods, 
in which the computational grid moves with the fluid, will not be covered at all. Time 
implicit methods will be mentioned, but not described in detail. Accelerated conver- 
gence to steady states and boundary stability theory will be ignored. Finally, the great 
majority of time explicit schemes, clever or otherwise, cannot be examined due to lack 
of time, especially as some of them are quite complicated. This is not a survey course in 
computational methods. 

What will be covered are the basic concepts fundamental to every CFD scheme. 
The emphasis will be on concepts and techniques which are simple, general, and have a 
clear mathematical basis (criteria which exclude a great many methods in use!). While 
the mathematical underpinnings are important, the usability of the methods is equally 
important; therefore mathematical rigor will not be attempted, as the focus will be on the 
concepts rather than their proof. Indeed, mathematically rigorous statements can seldom 
be made about solution methods for CFD problems. Virtually all rigorous mathematical 
work has been devoted to the simpler problems of linear systems or single nonlinear 
equations. Fortunately, techniques developed for these simpler problems usually work on 
the more complicated problems of nonlinear systems, even in the absence of mathematical 
proofs that they should do so. 

One standard of maturity for a field is the degree to which all conclusions follow 
logically from a small number of basic principles. By this standard, computational 
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fluid dynamics is not a mature field. It has been characterized by a large collection 
of complicated ad hoc methods. These lectures will try to achieve some maturity by 
getting a widely applicable approach from a small number of basic concepts. In this 
sense it departs from the “traditional” path in the CFD field, which has been to get the 
best possible solutions to specialized problems. 

Acknowledgment. I have performed this work while a National Research Council Re- 
search Associate, doing research into astrophysics in the Space Science Division at NASA 
Ames Research Center. I am indebted to the NRC, NASA Ames, and the Space Science 
Division for supporting me in this work. I would also like to express my appreciation to 
the NAS Project for sponsoring the production of the videotapes. 


Chapter 1 

Conservation Laws, Wave 
Equations, and Shocks 


This chapter deals with the general properties of hyperbolic systems of conservation 
laws, of which the most frequently encountered example is probably the fluid dynamics 
equations. It begins by defining what one means by “conservation,” and describes the 
close relationship between conservation laws and wave equations. We quickly discover 
that the concept of a continuous function as the solution of a partial differential equation 
is inadequate when we consider nonlinear equations, and are thus led to the concept of 
shock waves as the required discontinuous solutions. We will also find that not all shock 
waves which are allowed by conservation arguments are physically valid solutions, and 
that only those shocks which satisfy an entropy condition may actually occur. 


1.1 Conservation Laws 


The equations of fluid dynamics are one example of the mathematical formulation of 
conservation laws. The simplest one dimensional conservation law is a single equation, 
describing a single conserved quantity, and may be written as either an integral equation 
or as a partial differential equation. Suppose (in one dimension) that u is the density of 
some conserved quantity per unit length, and / is the flux of u (i.e., the rate at which the 
density u flows past a given x coordinate). The integrated density between two points 
Xi and xz , Xz < X\, satisfies [l] 


^ J u(x,t)dx = -\f{x u t) - f{x 2 ,t)}. 


( 1 . 1 ) 


If u and / are continuous functions of x and t, then in the limit as xz —■ ► xj = x, Eq. (1.1) 
becomes 


du df 
dt dx 


= 0 . 


( 1 . 2 ) 


Eq. (1.2) is said to be in conservation form. More generally, any equation in which the 
time derivative of a density plus the divergence of a flux equals inhomogeneous local 
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(i.e., non-derivative) source terms is in conservation form. In the case of (1.2), the flux 
divergence is df jdx, and the source terms are zero. 

The conservation form of an equation or system relates the time rate of change of 
a density in a small volume to the flux of that quantity through the boundary of that 
volume. Alternatively, one may recast the problem in terms of waves, and study the 
propagation of wave amplitudes. The two approaches are nearly the same when applied to 
single linear equations; however they are quite different, though still related, when applied 
to nonlinear systems such as the fluid dynamics equations. The equivalent wave equations 
are usually called characteristic equations in this context, while the wave velocities are 
called characteristic velocities. 


1.2 Wave Equations 

We begin by introducing the wave amplitude u, which is a function of the time t and the 
spatial coordinate x. Suppose there exists a curve C in the xt plane, parameterized by 
the variable o through the relations 


x = x{o), t = t(a), 


(1.3) 


where each point on C corresponds to a unique value of a. Then the rate of change of u 
along the curve C is [2] 


du du du 
da = St 1 ’ + HZ*’' 


(1.4) 


where the subscripts denote derivatives with respect to a. 
If u{x,t ) is constant along C, we have 



(1.5) 


and C is called the characteristic curve, or simply characteristic , of u. Eq. (1.5) is called 
a characteristic equation, and may be cast in different forms. 

A common form for the characteristic equation is obtained by dividing (1.4) (with 
du/da = 0) by t a (implicitly assuming that t a is never zero; we will assume that t{o ) is 
always an increasing function of a, as in t = a), giving 


du du 

al + a di =0 ' 


a — x a jt a — dxjdt , 


( 1 . 6 ) 


where a is the characteristic (or wave) velocity and is simply the slope of curve C in the 
xt plane. If a is constant, Eq. (1.6) is a linear wave equation. If a = a(u) constant, 
Eq. (1.6) is a nonlinear wave equation. 

The conservation and wave equations for u are equivalent provided that the wave 
speed a is given by 

_df 

a du' 


(1.7) 



1.3. CONTINUOUS ANALYTICAL SOLUTIONS 
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If a = constant, the problem is linear and the solution is 

u(x,t) = u 0 {x — at), where u(x,0) = u 0 (x). (1-8) 

The characteristics in the linear case are a family of parallel lines with constant slope 
dxjdt = a, along each of which u = constant, although u varies from one line to the 
next. The characteristics cover the entire xt plane in the region t > 0, implying that the 
solution exists as defined for all x and t > 0. 

The derivative du/do on the left side of (1.4) is often written with o = f, as in 

du du dx du ,, 

dt dt dt dx 

in which case it is called the total time derivative. The quantity dujdt in Eq. (1.9) 
represents the time rate of change of u as measured at a point moving with velocity 
dxjdt. 


1.3 Continuous Analytical Solutions 

Now consider the general case where a(u) is not constant. We know that u is constant 
along each characteristic curve C in the xt plane, hence a(u) is also constant and the 
characteristics are straight lines. However, as u varies from one characteristic to the 
next, so does a(u), and the slopes of the characteristics also vary from one to the next. 
Those characteristics which are initially separate will either diverge or converge, leading 
to behavior not seen in the linear problem. A typical characteristic diagram is shown in 
Figure 1,1, reproduced from page 21 of the book by Whitham [l]. 

Let the initial value problem be given by Eq. (1.6) with u(x,0) = /(x). Then for 
each characteristic curve which intersects the x axis at x ~ £ at time t = 0, u = /(£) 
everywhere on that curve. The slope of a given curve is a(/(£)) = F(£) and is also known. 
Each characteristic is therefore defined by a unique £ value, and has as its equation 

x = Z + F(Z)t. (1.10) 

The complete set of £ values defines a whole family of characteristics, along each of which 
the solution is 


II 

s 

II 

3 

II 

«(/(£))■ 

(in) 

Now let’s check the solution. We find 




1 

d£ 

f(e) 

(1.12) 

dx 1 + F'{t)t’ 

dt 

1 + F'{t)t' 


where the primes denote derivatives with respect to £, and therefore 

du _ E(fl/'(fl du _ /'(Q 

dt l + F'{£)t ’ dx 1 + F'{£)t ’ 

and Eq. (1.6) is satisfied. 


6 


CHAPTER 1. CONSERVATION LAWS, WAVE EQUATIONS, AND SHOCKS 


CONTINUOUS SOLUTIONS 



Figure 1.1: Characteristic Diagram for Nonlinear Waves. Copyright © 1974 by John 
Wiley & Sons, Inc. Reprinted by permission of John Wiley & Sons, Inc. 

Saying that u is constant along a characteristic curve whose slope is a(u) = dxjdt 
is equivalent to saying that each particular value of u propagates at a characteristic 
velocity a(u), which is the wavelike behavior we were seeking. In the linear case where 
o(u) = constant, the solution obtained above reduces to the simple form of Eq. (1.8). 
In the nonlinear case, the velocity changes from point to point and the wave exhibits a 
nonlinear distortion. 


1.4 The Breaking of Waves 

Difficulties arise when the wave undergoes compression, which occurs in any region where 
the velocity a is decreasing, i.e., in any region for which F'(£) < 0. Any two character- 
istics in this region which are initially separate will cross at some later time, giving a 
solution which is multivalued. This difficulty is apparent in the solution of (1.13): the 
wave “breaks” (acquires an infinite slope) at the time t = -1/F'{£). Breaking occurs 
earliest on the characteristic defined by £ = for which F‘(£) < 0 and |F'(£)| is the 
maximum, and at a time given by 


tB f'UbY 


(1.14) 


An extreme case of breaking occurs when the initial state is a step function at x = 0, 
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tL\m d\ — ^ / -g -i e\ 

u(z,0) = f(x) = < a(x,0) = F(x) — < (1.15) 

u 2 , a 2 = a(u 2 ), x < 0. 

If a 2 > ai the solution breaks immediately. 

A different result occurs of a 2 < a x . The wave undergoes expansion, and a continuous 
solution results. At times t > 0 all values of -F(O) between a\ and a 2 spread out along a 
fan of characteristics passing through £ = 0. This expansion fan is continuous, and must 
satisfy (1.10) and (l.ll), hence 


X X 

a — a 2 < — < ai. 


(1.16) 


and the complete solution for a is 


a j , dj ^ x j t , 
a = \ x/t, 0-2 < x ft < Oi; 

flj , X j 1 ^ fl 2 . 


(1.17) 


The solution for u(x,t ) is then obtained by inverting the known relationship a — a(u). 


1.5 Shock Waves 

In practice, the nonlinear wave equation we wish to solve represents a process whose 
physical reality must be single valued, and the multi-valued solution produced by wave 
breaking must be rejected. Yet solution (1.10)— (1.13) is valid up until the derivatives 
become infinite, so we must modify our concept of a solution to include discontinuous 
solutions which are single valued, and have a finite number of discontinuities. The formu- 
lation of the wave problem as a partial differential equation is not valid for discontinuous 
solutions, because the derivatives are defined only for continuous functions. However, 
the integral formulation of conservation law (1.1) is valid even when u is discontinuous, 
and it is to this form we turn now. 

Stable discontinuities in nonlinear waves are called shock waves, or shocks. On either 
side of the shock the solution is continuous and differentiable. If the shock is located at 
position x,{t) and moves with velocity v t , where x 2 < x„(t) < xi at time t, we may write 
the integral equation for u as [l] 

/(*a»*) - /(*i*0 = ^J z u(x,t)dx+ J u(x,t)dx, (1.18) 

/ x7 Q rxi Q 

— u(x,t) dx + u{x~,t)v t + J + —u(x,t) dx - ti(x+,t)tv(1.19) 
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In the limit as x 2 — » x t , X\ — » ► x+ this becomes the shock jump condition 

fi ~ fi = v t (u 2 - ux), (1-20) 

where /i, u\ are the flux and value of u on side 1 of the shock (taken to be the right 
here), and / 2 , u 2 are the values on side 2. 

The shock jump condition may also be obtained by transforming the conservation 
equation (1.2) to a reference frame in which the shock is at rest; i.e., which moves at a 
velocity v 8 with respect to the original coordinate system. Integrating the transformed 
equation over a small volume around the shock, and taking the limit as the volume goes 
to zero, shows that the transformed flux / — uv s must be continuous across the shock. 
Thus /i — Ui v 9 = f 2 — which gives Eq. (1.20). 

Rewriting (1.20) gives the shock velocity as 

fM - /(« 0 

V 8 = — 

U 2 — u 1 

In the weak shock limit where |u 2 — t*i| |u 2 | + \u\\ we have u 2 « U\ = u, and v 8 — ► 
d f fdu = a(u). Thus a sufficiently weak shock is a small discontinuity which travels at the 
local wave velocity. A strong shock, however, has a velocity which is distinct from both 
u\ and u 2 . Note that v 8 = a always for the linear wave equation, where a = constant. 



1.6 The Inviscid Burger’s Equation 


Perhaps the most well known nonlinear wave equation is the inviscid Burger’s equation 

du 


dt 


+ 


du du d f u 2 \ 

u ^ = 0 ’ or sr + a;(yj =0 ' 


( 1 . 22 ) 


which has f(u) = u 2 / 2, a(u) = u. Burger’s equation appears in most texts on nonlinear 
hyperbolic equations, and in the inviscid form represents the velocity field of a gas of 
non-interacting particles (such as a cloud of dust particles in a vacuum, or a zero pressure 
gas). 

The inviscid Burger’s equation satisfies the general solutions obtained above. Its 
solution for an initial discontinuity at x = 0, with u = u 2 at x < 0, u — ui at x > 0, is 
an expansion fan for u 2 < Uj: 


u 2 , 

xft < u 2 \ 

x/t, 

u 2 < x/t < «i; 

U 1 , 

ui < x/t ; 


( 1 . 23 ) 


and a shock wave for u 2 > u 2 : 


1.7. EXPANSION SHOCKS AND THE ENTROPY CONDITION 
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where the shock velocity is obtained from (1.21): 


», = 1(U, + u,). (1.25) 

The linear wave equation and the nonlinear Burger’s equation represent special cases 
of one dimensional planar fluid flow. If the fluid velocity v is constant with x, then 
the continuity equation for the density reduces to (1-6), with u the density and a = v = 
constant the fluid velocity. Density inhomogeneities are carried along by the flow without 
distortion. However, if the pressure is spatially constant ( dp/dx = 0), then the velocity 
equation reduces to the inviscid Burger’s equation, with u the fluid velocity. 


1.7 Expansion Shocks and the Entropy Condition 

The preceding analysis found two different solutions to the nonlinear wave problem with 
a step function initial condition. The case with a 2 > ai remains a step function, the dis- 
continuity being a shock wave which propagates at the shock velocity given by Eq. (1.21). 
The case with a 2 < a\ breaks up into an expansion fan (which widens with time) adjoining 
the constant states u 2 and u i on either side. 

The shock jump condition (1.20) says nothing about the breakup of a discontinuity 
into an expansion fan. It is natural to ask if a discontinuity between two states with 
a 2 < a\ could also propagate as a shock, since the possibility is permitted by the jump 
condition. Such a shock is known as an expansion shock , since it leaves behind a state of 
decreased “density” u. 

It turns out that expansion shocks are unstable. Small perturbations to an expansion 
shock solution grow, destroying the shock. Compression shocks, on the other hand, are 
stable. Stable shocks must satisfy the entropy condition [l], [3] 

a 2 > v s > ai, (l- 26 ) 

¥9 i * • Cl 

which means that characteristics cross at the shock front, and the slope of the shock 
front line in the xt plane lies between the slopes of the characteristics which intersect 
it from either side. Consequently no characteristic drawn in the direction of decreasing 
t intersects a line of discontinuity, and every point in the plane can be connected by a 
backward drawn characteristic to a point on the initial line at t = 0. In the case of the 
inviscid Burger’s equation, a = u, and v 3 = (u 2 + Ui)/2, which satisfies (1.26) only for 
u 2 > Ui, as expected. 

Expansion shocks are forbidden solutions for systems of nonlinear equations as well. 
In the case of fluid dynamics, an expansion shock would cause the entropy of the flow to 
decrease with time, and is forbidden by the laws of thermodynamics as well as stability 
considerations. It is for this reason that criteria such as (1.26) are known as entropy 
conditions. 
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1.8 Nonlinear Systems in One Dimension 


In this section we consider systems of conservation laws 


oi me iorm 


du df 
dt dx 


= 0 , 


u = 


A 



u l 


r 

f 1 

u 2 


f 2 


, f = 


U " J 


< fn J 


(1.27) 


where f is a vector flux and a function of the conserved densities which are the components 
of the vector u. (Note that the term vector is used here in the linear algebra sense to 
represent any set of unknowns, and not, for example, as a velocity vector.) 

As in the previous section, let the states on either side of a shock be numbered 1 and 
2 (1 for the right side, 2 for the left). The shock jump conditions are obtained by writing 
Eq. (1.27) in integral form, and integrating over a small region containing the shock, as 
in section 1.5. The n jump conditions are therefore 


fa ~ fn = v,(u i2 - u.i), * = 1 , . . . , n, (1.28) 

where fn — /i(ui) is the ith flux on side 1 of the shock, and similarly for /,- 2 . The shock 
velocity v t is the same for all n jump conditions. 

Eq. (1.27) can also be written in the form 


du du dui ,A duj 

+ A— -°, or — + - °, t-l,...,n, 


dt 


(1.29) 


where A is the Jacobian matrix of elements a,j defined by 


dL 

a„ = (1.30) 

The system is defined to be hyperbolic if the matrix A has n real eigenvalues A,, 
which we will order so that Aj < A 2 < • • ■ < A n . The eigenvalues are the characteristic 
velocities at which signals propagate in the medium described by (1.27). The character- 
istic. velocities are the slopes of the characteristic curves, as was the case for the single 
wave equation. However, constancy of any wave amplitude along a characteristic curve 
no longer implies constancy of the wave velocity along that curve, so the characteristic 
curves need no longer be straight, and in general will not be. (A more detailed discussion 
of the characteristic equations for nonlinear systems will be given in Chapter 5.) 

The characteristic velocities determine an entropy condition for the nonlinear system, 
which defines the class of allowed shock waves in a fashion analogous to Eq. (1.26). The 
entropy condition for the nonlinear system, as given by Lax [3], is that for some value of 
k, 1 < k < n, the inequality 


A*(u 2 ) > v, > A*(ui) 


(1.31) 


1.9. EXERCISES 


11 


must be satisfied. 

In the case of fluid dynamics, the characteristic velocities are Ai = v — c, A 2 = v, and 
A 3 = v + c, where v is the fluid velocity and c is the speed of sound. If the shock speed 
is positive, Eq. (1.31) is satisfied by the case k = 3: 

v 2 + c 2 > v, > V! + Ci, (v, >0), (1.32) 

which means that right-moving sound waves on either side of the shock intersect the 
shock. If the shock speed is negative, then 


v 2 -c 2 >v s > t>! -ci, (u, < 0), (1.33) 

which is the k = 1 condition, and implies that all left-moving sound waves intersect the 
shock. 

In either case, sound waves behind the shock (in the compressed region) travel faster 
than the shock and catch up to it, while sound waves in front of the shock are propagating 
more slowly than the shock and are overtaken by it. Thus the shock velocity is subsonic 
relative to the post-shock state, but supersonic relative to the pre-shock state. Very 
weak shock waves (characterized by very small shock jumps) are simply sound waves, 
and travel at the speed of sound relative to the fluid. In this limit, the three velocities 
in (1.32) become the same. 


1.9 Exercises 

1. Verify that the solution of (l.lO)-(l.ll) satisfies the integral equation (l.l), for the 
linear case where F(£) = a = constant. 


2. Find the shock jump condition for the equation 


dt 




= 0 . 


Is it the same as (1.25) for Burger’s equation? Does this result seem unusual? 

Comment: A particular conservation law, such as (1.22), will give rise to an infinite 
number of other equations upon multiplication by u to some power. The resulting 
equations are equivalent for u a continuous function, but give rise to different shock 
jump conditions. The correct jump condition may only be determined by reference 
to the physics of the problem, i.e., by working with those quantities which are both 
conserved and physically meaningful, such as mass, momentum, and energy. 


3. The general expansion fan u = (x — Xo )/(t — to) is a solution to the inviscid Burger’s 
equation (1.22). Consider an initial value problem given by u = x/(l + t) for x < xo , 
u = 0 for x > xo, at time t = 0. The point xo is the initial shock location. Find the 
solution for t > 0, including the shock location x, and shock velocity v,. (Hint: v, 
is not constant with time.) 
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4. The fluid equations for a perfect gas in one dimension are 


dp 

dt 


+ Tx {pv) = 0 , 


d d . , 

alW + 


d 

(\ , 'l 

d r 

/1 , \ 1 

dt 

li'”’ + V 

l H 
+ 

l-pu + e + pi v 


where p is the mass density, v is the velocity, p = (7 — l)e is the pressure, e is the 
thermal energy density, and 7 is the (constant) ratio of specific heats (= 5/3 for 
an atomic gas). Find the shock jump conditions. Show that 


Pi _ (7 + l)p 2 /px + (7 - 1) 
Pi (7 - IWpi + (7 + 1 )’ 


If P 2 / Pi can range from one to infinity, what values can the density ratio pij pi take 
on? (Hint: Although there are now three shock jump conditions, there is still only 
one shock velocity.) 

5. Using the shock jump conditions derived in the previous problem, let Pi = pi = 1 , 
v i = 0> 7 — 5/3, and p-ijpi = 5. The speed of sound is c = yj'ip/p. Compute p 2 , v 2 , 
and v ,, draw a diagram of the pre- and post-shock characteristics, and show that 
condition (1.32) is satisfied. 

Condition (1.32) is condition (1.31) with k = 3. For this problem, is condition 
(1.31) satisfied for any other k? 



Chapter 2 

Finite Difference Solutions of Wave 
Equations 


The first chapter considered the general properties of conservation laws and wave equa- 
tions in one dimension. This chapter will cover the basic concepts involved in the approx- 
imate solution of wave equations by finite difference methods. The conservation laws of 
interest are partial differential equations in one time and one or more spatial coordinates. 

The numerical techniques to be described in this chapter will be oriented at first 
toward the solution of the linear wave equation, in order to introduce the basic concepts 
without obscuring them by the complexity encountered in solving nonlinear systems. 
However, since the solution of nonlinear systems in more than one spatial dimension is 
our ultimate goal, the implications of nonlinear systems for the solution techniques pre- 
sented will be discussed after each technique is described. In Chapter 4 we will find that 
consideration of the very general problem of multidimensional nonlinear systems with 
shocks will determine the form of an artificial viscosity. This chapter will demonstrate 
the need for artificial viscosity in problems with discontinuous solutions and will intro- 
duce a form suited for the problems under discussion, but the fundamental justification 
of the form will be deferred to Chapter 4. 


2.1 Basic Finite Difference Approximations 

We begin by considering a piecewise continuous function f(x ), defined for all x. Ideally, 
we would like to know the exact value of / for all x values, but in practice we will be 
restricted to knowing / at a finite number of discrete points. Describing a function in 
terms of its values at discrete points is known as discretization. For convenience sake, we 
will discretize f(x) by recording its values at the set of grid points x,-, spaced a uniform 
distance Ax apart, 

x, = t'Ax, (2-1) 

although in general one can have arbitrary grid spacing, which leads to more complicated 
finite difference approximations. Denote the values of / at Xj by /<: 

U = /(*). ( 2 - 2 ) 


13 


14 


CHAPTER 2. FINITE DIFFERENCE SOLUTIONS OF WAVE EQUATIONS 


Given the fi values, we need to approximate various derivatives of / at the grid points. 
The starting point for all such approximations is the discrete Taylor’s series [4], 


fi+n = fi + n Ax 


df_ 

dx 


+ 


(nAx) 2 d 2 f 
2 ! ~dtf 


+ ••■ + 


(nAx) m d m f 
ml dx m 


+ 


(2.3) 


which holds for /(x) a continuous and infinitely differentiable function. If /(x) is not 
continuous, this approximation breaks down and difficulties may occur. Techniques for 
dealing with discontinuous functions will be discussed later. For now we will assume that 
all functions are continuous. 

Setting n = ±1, we see immediately that 


df_ 

dx i 


^ (/«+ i - fi) + O(Ax), 
Ax). 


(2.4) 

(2.5) 


Eqs. (2.4) and (2.5) are called one sided difference approximations, and are the basis 
for the low order upwind schemes to be discussed later. These one sided formulas are 
formally of first order accuracy. The order of an approximation is given by the power 
of the grid spacing Ax appearing in the leading error term, since the error in the ap- 
proximation vanishes as that power of Ax in the limit as Ax — »• 0. We can see that the 
one sided approximations are first order accurate by solving for df /dx| f - in Eq. (2.3) with 


n = 1: 

df_ 

dx 


1 

Ax 


(/«'+! ~ fi) ~ 


Ax d 2 f 


(Ax) 2 d 3 f 
3! dx 3 , 


(Ax)™" 1 d m f 
m! dx m i 


( 2 . 6 ) 


For a sufficiently small Ax, the leading error term dominates, hence Eq. (2.4) is a first 
order approximation to df/dx\i, as is Eq. (2.5). 

The first order approximations are useful in certain situations (such as at a bound- 
ary, where symmetrically located data are not available), but their slow convergence 
properties make them undesirable in most instances if more accurate approximations are 
available. 

Replacing n by ±n in Eq. (2.3) and taking the difference of the expressions yields the 
following result: 


■2 fi+n 


rp2 — 
1 n — 


U 


2nAx 


df_ 

dx 


+ 


(nAx) d f 


3! dx 3 


(nax J’ cr/ 


5! dx 5 


+ 


U J 


7! dx 7 


+ (2.7) 


from which we get the second order approximation 


d[_ 
dx i 


/«+! ~ fi-1 
2Ax 


+ 0(Ax 2 ). 


( 2 . 8 ) 


Note that Eq. (2.8) is not the only centered second order approximation. The quantity 
T 2 , defined by Eq. (2.7), is a second order approximation to d//dx|, for any n value. 
However, the error is quadratic in n, so the case n = 1 is clearly the best choice. 
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We’ve seen first and second order approximations to df jdx. Higher order approxi- 
mations may be constructed by taking appropriate linear combinations of the formulas 
already obtained. For example, Eq. (2.7) gives an exact expression for the error terms 
in the approximation T 7 . The following linear combination gives a fourth order approx- 
imation to df jdx\i". 


4 - ~ T n 


T* = 


n 4 - n 2 Az 4 d 5 f 


ri 6 — n 2 Az 6 d 7 f 


n 2 — 1 dx\i n 2 — 1 5! dx h I, n 2 — 1 7! dx 7 1« ^ 

The leading error term is again quadratic in n, and is minimized for n = 2 (note that n 
may not be 1), giving 

T} = jT? - \t J (2.10) 

as the best fourth order approximation to df /dz|,-: 


df_ 

dx 


[8 (/,•+ 1 — U- 1 ) — (/.+ 2 — /*— 2)3 + 0(Az 4 ). 


( 2 . 11 ) 


12Az 

Higher order approximations are obtained by computing lower order approximations 
over different intervals nAz, then taking linear combinations of the results with coeffi- 
cients chosen to cancel out the leading error term. 

The above approach is an application of Richardson extrapolation , or the deferred 
approach to limit [4]. Richardson extrapolation can be used to obtain improved, higher 
order estimates of any quantity whose error is known to consist of a power series in some 
discretization parameter. One computes the lowest order approximation to the quantity 
using different discretization parameters (such as nAz in the above example), then takes 
linear combinations of the two most accurate lower order approximations to obtain a 
higher order approximation. In the above example, only even powers of the discretiza- 
tion parameter appeared in the power series, so that a fourth order approximation was 
obtained very quickly. In other circumstances, such as the first order approximations of 
Eq. (2.6), all powers of the discretization parameter appear and more work is required 
to obtain a given order of accuracy. 

Another common example of Richardson extrapolation is Romberg’s method for ap- 
proximating definite integrals. One computes a sequence of trapezoidal approximations 
to the integral, using different interval sizes, then takes linear combinations of the re- 
sults to eliminate the leading error terms. The resulting approximation yields a much 
more accurate approximation for a given amount of effort than does the trapezoidal 
approximation by itself. 

A similar analysis leads to the following formula for the second derivative of a function, 
to second order: 


d'f 


dx 2 


Ax' 


:(/»+! - 2/,- + ft- 1 ) + O (Az 1 ). 


One also encounters diffusion terms, which may be approximated as 


Uf\ = /.(/(« - fi) ~ «,-./*(/< - /.-.)] + 0( Ax J ), 

and which reduces to Eq. (2.12) for the case k = 1. 


A 

dx 


( 2 . 12 ) 


(2.13) 
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2.2 A Brief Survey of “Traditional” Solution Meth- 
ods 

Having mastered the basics of finite difference approximations, we will now attempt to 
solve the one dimensional linear wave equation 

du du 

— + a— — 0, a = constant, (2.14) 

on a grid in the xt plane, at the points (x,-, t n ), spaced evenly in the x direction (x,- = t'Ax) 
but' not necessarily in the t direction. Let U? = U(xi,t n ) be the numerical approximation 
to the exact solution u(x,i ) at the grid points. 

We are given a partial differential equation which relates the time variations in a 
quantity u to its spatial variations, as well as initial data for all x values. The approximate 
solution at a time At later is obtained from a numerical approximation to the differential 
equation. Repeated applications of the numerical approximation yield solutions at a 
sequence of advancing time steps. 

The solution to the problem is completely defined once we have the differential equa- 
tion and the initial data; hence the problem is called an initial value problem. A related, 
but more complicated, problem is the initial boundary value problem , in which boundary 
conditions as well as initial conditions over a finite spatial domain are specified. I will 
return to the issue of boundary conditions in Chapter 5, but for now we will look at the 
finite difference solution of wave equations over an infinite domain. 

The literature describing methods for solving wave equations and conservation laws is 
very large, and a comprehensive examination of all such methods is not feasible. Instead 
we will look a few representative and commonly used methods, chosen not only because 
they illustrate the basic concepts of numerical solutions but also because they may be 
applied to a wide variety of problems. 


2.2.1 First Order Upwind Methods 


The simplest and least accurate solution methods for the linear equation are the first 
orders “upwind” methods. Replacing the time and space derivatives in Eq. (2.14) by their 
first order approximations yields four possible solutions: 


where 


U\ 


u ? +1 = 

17" +1 = 
r " +1 = l 

n+1 j 


v? - o{U? - l/JLj), 

(2.15) 

u?-o(u? +l -un, 

(2.16) 

?-o{ur i -v&). 

(2.17) 

? -"{Vm - V? +1 )> 

(2.18) 

At 

(2.19) 

a = — 

Ax 


is called the Courant number , after Richard Courant, whose work in linear and nonlinear 
waves and their solution forms the foundation for much of the field today. 
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At first glance these four approximations may seem equally valid-and they are from 
the standpoint of the formal accuracy of the approximations out of which they were 
constructed. However, two of them are guaranteed to be unstable for any time step 
size At . Stability analysis will be deferred to Chapter 3, but for now it will simply be 
stated that Eqs. (2.16) and (2.18) are unstable if a > 0, while Eqs. (2.15) and (2.17) are 
unstable if a < 0. “Instability” means that the numerical solution grows exponentially 
with the number of time steps, even though exponential growth is not a valid solution 
to the given initial value problem. 

The stability criteria for the upwind methods illustrate a very general property of 
numerical solution methods, which can be stated as follows: The domain of dependence 
of the numerical approximation to the solution of the differential equation must include 
the domain of dependence of the original differential equation . For example, a right- 
moving wave (with a > 0) has a solution whose time evolution at a particular point 
is governed by the spatial variation of the wave to the left of that point— not to the 
right. Therefore the numerical solution must make use of information to the left of the 
grid point, which Eqs. (2.16) and (2.18) do not. One can define more precisely which 
grid points should contribute to the solution by examining the characteristic curves of 
the differential equation, but for now it will be sufficient to state that one may include 
more points than are actually required, perhaps to improve accuracy, but the stability 
requirements will depend in detail on the particular scheme chosen. Note that the upwind 
scheme is so named because the points to be included in the first order approximation 
are “upwind” from the current grid point. 

Assuming that a > 0, we find that Eq. (2.15) is stable provided 0 < a < 1, while 
Eq. (2.17) is stable for any u > 0. If the largest stable \a\ is a maz , we may use any At 
value satisfying 

( 2 . 20 ) 

N 

Thus approximation (2.15) is stable for At < Ax/|a|, while (2.17) is stable for any A£. 
Similar stability criteria apply to approximations (2.16) and (2.17) when a < 0. 

2.2.2 Explicit vs. Implicit Methods 

The approximation of Eq. (2.15) is an example of a time explicit , or simply explicit , 
method. An explicit method is one in which an unknown value at time step n + 1 appears 
at only one grid point in the formula (e.g., J7^ +1 ), and is given explicitly in terms of the 
previous values at step n. The approximation of Eq. (2.17) is an example of a time 
implicit , or simply implicit , method. An implicit method is one in which the unknown 
values of U at step n + 1 appear at more than one grid point in the approximation, and 
are determined implicitly by a set of simultaneous equations, rather than by a set of 
independent explicit equations. 

Explicit and implicit methods each have their advantages. Explicit methods are 
much simpler to implement, as they do not entail the solution of simultaneous equations. 
Explicit methods usually require much less time to evaluate, per time step, than implicit 
methods, because the simultaneous equations inherent in the implicit schemes usually 
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require lengthy calculations to solve. On the other hand, implicit methods are stable 
for much larger time steps (often infinite, as above) than explicit methods. However, 
the error in the time dependent numerical solution increases as some power of the time 
step, so that in general an implicit scheme would require roughly the same time step as a 
similar explicit scheme to achieve the same level of accuracy, but at a much higher cost in 
computer time. Thus explicit methods are almost always preferable for time dependent 
problems. 

Nevertheless, there are situations in which implicit methods are preferable. One case 
is the solution of “stiff” problems, which contain large characteristic velocities whose 
corresponding effects are of minor importance. The time step for an explicit scheme is 
limited by the largest characteristic velocity, whether or not the associated phenomena 
play a major role in the solution. An implicit scheme may then be used with a time step 
small enough to follow the phenomena of interest (such as convection), but much larger 
than that dictated by the time scales of unimportant phenomena (such as sound waves). 
Flows at very small Mach numbers fit into this category. 

Another case in which implicit methods are useful is in steady state flow models. One 
way to model a steady state flow problem is to pick a set of initial conditions and let the 
problem evolve to an equilibrium state. The time dependent solution is of no interest, 
and need not be accurate provided that the steady state solution obtained is correct. 
Once again, a small number of computationally expensive steps may be more efficient 
than a large number of inexpensive steps. 

Implicit methods generally have to be tailored quite carefully to the problem at hand, 
while explicit methods can be made very general. Thus from now on the discussion will 
focus on explicit methods. 


2.2.3 Two Popular Second Order Schemes 


The first order upwind schemes suffer from two deficiencies. The first is low resolution. 
An initially sharp profile (such as a step function) is smeared out or “diffused” over many 
grid intervals as the solution progresses. This smearing tendency is known as numeri- 
cal viscosity, and is a nonphysical effect (nonphysical because the original differential 
equation was inviscid). 

The second deficiency concerns the nature of the upwind stability criteria. In subsonic 
fluid flow there is no simple upwind direction, as the flow has characteristic velocities 
in all directions. A direct application of the upwind scheme will therefore be unstable, 
unless the equations are put in characteristic form first, and each characteristic equation 
approximated by the appropriate upwind form. The characteristic form is useful for 
specifying boundary conditions, but is needlessly complex for other problems. 

Finite difference schemes which make use of symmetrically located data points possess 
stability criteria which are independent of the flow direction. One early attempt at a 
symmetric method, which is still widely used, is the Lax-Wendroff scheme [5]. We start 
with a Taylor’s series approximation for U* +1 about time t n , truncated after the second 
order term: 



At 2 d 2 U 
2 dt 2 


n 

i 


= U? + At — 


( 2 . 21 ) 
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Equation (2.14) may be solved by substituting — df jdx for dU/dt in (2.21), and using 
Eq. (2.5) for the wave velocity. The result is the approximation 


u; 


1 __ 


U? - A t 


df_ 

dx 


+ 


At 2 d 
2 dx 



( 2 . 22 ) 


where a = a(U) in general. Making the spatially centered, second order finite difference 
approximations 


d£ 

dx 


2Ax 


(fi+i- fi-x). 


(2.23) 


d_ 

dx 



Ax 2 


[ a » + l/2 (/»'+! — ft) a i-l/i{fi > 


(2.24) 


gives an explicit solution to Eq. (2.14) which is second order accurate in time and space. 
When applied to the linear wave equation, the Lax-Wendroff method gives 


t/“ +1 =U?- i<r(lC, - U ?_ ,) + - W? + £?_,), 


(2.25) 


with a given by Eq. (2.19), which has a stability criterion independent of the sign of a: 

< 1. (2.26) 


At 

a-.~ 

Ax 


The Lax-Wendroff method is easily applied to single one dimensional wave equations. 
However, its extension to systems in one dimension, with an unknown variable vector 
u and vector flux function f(u), entails the calculation of the Jacobian matrix di j dM 
where the innocent-looking a appears in Eq. (2.22), a tedious if not impossible task. 

MacCormack [6] came up with an apparently simpler method for the nonlinear sys- 
tems case. MacCormack’s method reduces to Eq. (2.25) for the linear wave equation, 
but is not the same as the Lax-Wendroff method for more complicated problems. The 
MacCormack method solution to Eq. (2.14) may be written 

; . ;■ ;\i l «: : 

u< = !7”-|^(/r +1 -/n. 

v r 1 = u? - - rn, (*•») 


1 

2 


(or + ft) - 


At 

2Ax 


Cf> - /-.)• 


This is a two step explicit method which does not require the Jacobian of the flux 
function. However, it is not a symmetric method: Eq. (2.27) uses a right direction one 
sided approximation to df n jdx and a left direction one sided approximation to df jdx. 
One could have chosen the left direction approximation for the derivative of f n and 
the right direction for f. The choice of direction for the terms is arbitrary, so long as 
they are opposed, but the two possible choices will give slightly different results when 
applied to the same nonlinear problem. MacCormack advocated switching the direction 
of differencing at successive time steps in order to restore symmetry. Switching entails 
cycling through 2 m possible schemes in succession in m dimensions, which can lead to 
complications as difficult as those of the Lax-Wendroff scheme. 
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The Method of Lines 


2.3 


The preceding methods, while simple to formulate for single one dimensional wave equa- 
tions, become quite complex and tricky to implement for multidimensional nonlinear 
systems such as the fluid equations. The difficulties stem from the use of one-sided 
difference approximations in the case of the upwind schemes, and from the simultane- 
ous approximations of space and time derivatives in the Lax-Wendroff and MacCormack 
methods. These difficulties can be avoided by the method of lines approach, which is 
readily applicable to any time dependent partial differential equation solution. 

We are given a set of partial differential equations in space and time, along with the 
initial data at time t = 0. Approximating the spatial derivatives with finite difference 
®^prcssions in turn tells us how the solution changes in time at each grid point, allowing 
us to integrate the time derivatives to obtain the solution at a new time step. 

More formally, suppose we have a system of equations describing the components of a 
solution vector u. (For example, we might have u = ( p , m, e), where p is the fluid density, 
m the momentum density, and e is the energy density.) If the time derivatives of the 
components Uj appear only in the first degree, then the system can always be written as 


du 

~dt 


= P u, 


(2.28) 


where P is an operator involving any combination of the coordinates x and t, as well as 
the unknown variables u, and their spatial derivatives, but no time derivatives of u ; . In 
the case of Eq. (2.14), u is a vector with only one component u, and Pu = —df / dx. 

Now approximate all spatial derivatives in .Pii with the appropriate finite difference 
formulas, yielding the numerical approximation (PU),- (e.g., a spatially centered approx- 
imation such as (PU)i = -(/ l+1 - fi-i)/2Ax). Then we have a semi-discrete equation 
for dVi/dt, the time derivative of U at grid point t: 


~df = {pv) ‘- < 2 - 29 > 

In principle, any ordinary differential equation solution technique may be used to solve 
the coupled set of ordinary differential equations in Eq. (2.29), as long as the stability 
properties of the algorithm chosen allow reasonable step sizes At. Many such methods 
are available. 

The reduction of a system of partial differential equations to semi-discrete form, 
followed by the integration of the system by an ordinary differential equation solver, is 
known as the method of lines. 

One effective integration method for Eq. (2.29) is the classical, four step, fourth order 
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Runge-Kutta method [4], written as: 

Uj 1} = U?+ |A<(PU n ),-, 

ui 2 > = v T +\Atm< : {2 30) 

uj 3) = U? + At(PU^j<, 

U? +1 = U? + iA<(PU n + 2PUW + 2PU( 2 ) + PU( 3 )),-. 

This method requires the storage of the ofd step tJ”, the current intermediate step , 
and the running total of the sum of operators appearing in the last step. 

The time error in (2.30) is 0(At 4 ). Since the maximum At is proportional to Ax 
because of the stability criterion (given below), the time integration contributes an error 
of 0(Ax 4 ). Thus the spatial approximations should also be fourth order accurate, or 
else some of the effort involved in the time integration is being wasted. If the spatial 
approximations are at best second order accurate (common in fluid problems), we should 
use a second order time stepping scheme which requires less work to perform than (2.30). 
The following four step, second order method is suited for such problems [7]: 


Uj x) = U? + \At{PV n ) it 

U,- 2) = U? + §Af(PUU)i, 

U[ 3) = U? + f Af(PU< 2 )),-, 

U? + * = ' U? + At(PU (s) )i. 


(2.31) 


Not only does Eq. (2.31) require less storage than Eq. (2.30), but its simpler structure 
allows all four steps to be performed by one master loop, using a different coefficient for 
At in each iteration. 

Both the fourth order Runge-Kutta method of (2.30) and the second order method 
of (2.31) have the stability criterion 


\cr\ = 

At 
a — — 


Ax 


< o. 


_ 2 : w max 


(2.32) 


for the linear wave problem. If a second order centered approximation is made to du/dx, 
then cr max = 2\/2; if a fourth order centered approximation is made, then o max = 2.06. 
(Note: when Pu is a linear operator, such as a du/dx, with a = constant-and only 
then-Eq. (2.31) also gives a fourth order accurate time integration.) 


2.4 Dissipation and Discontinuous Solutions 


Consider the initial value problem 

~ + a— = 0, a = constant > 0; (2.33) 

dt dx 
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. 1 x < 0, 

u(x,0) = u 0 {x) = < (2.34) 

[ 0 x > 0. 

The solution is a step function of unit height, moving with velocity a. We would like to 
solve this problem using the finite difference methods discussed so far. 

This innocuous-seeming problem is in reality one of the most difficult wave problems 
to solve numerically. The techniques previously described give uniformly poor results at, 
say, a Courant number of o = 0.5. The results fall into one of two categories: 1) solutions 
which are monotonic but badly smeared, with the original discontinuity spread out over 
many grid intervals, and 2) solutions which contain spurious oscillations, especially near 
the discontinuity, but where the discontinuity is more sharply resolved than in category 
1. The first order (upwind) methods yield the diffuse monotone solutions, while the 
higher order methods yield the oscillatory solutions. 

These results are explained by the following theorem [8]: Any linear finite difference 
scheme which is guaranteed to preserve monotonicity is no more than first order accurate. 

This theorem is a major disappointment, as we would like to create high order solution 
methods which preserve the monotonicity of discontinuous solutions. The diffuse solu- 
tions produced by first order methods are usually inadequate; however, highly oscillatory 
solutions for discontinuous problems are equally inadequate. The theorem states that no 
linear finite difference scheme will satisfy the conflicting requirements of monotonicity 
and resolution. Therefore we must turn to nonlinear schemes for improvement. While 
the properties of nonlinear schemes can seldom be analyzed analytically, much can be 
accomplished by the careful blending of experience and linear theory. 

We begin by examining more closely the first and second order approximations in 
semi-discrete form: 

dUi 1 

(First order), (2.35) 

= “ a* 2^‘ +1 ~~ ( Secon< * order )- (2.36) 

The choice of time integration is irrelevant for the moment, but Eq. (2.31) will do in 
practice. Note that Eq. (2.35) may be written as Eq. (2.36) plus a correction term: 

dU{ 1 o' I o' I 

^ ~ 2 ^ Ui + l ~ u i- 1 ) + g (^ +1 - 2Ui + U'- 1 ) (First order). (2.37) 

(The above equation is in fact a general upwind scheme, which automatically switches to 
the appropriate direction based on the sign of a.) The correction term may be thought 
of as the finite difference approximation to the dissipative term 

d ( |oj Ax 2 du\ 
dx \ 2 At dx) ’ 

with a diffusion coefficient defined by the Courant number o and the grid parameters 
Ax and At. The added dissipation damps (i.e., prevents) the spurious oscillations which 
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would otherwise be generated by the second order centered difference approximation to 
du/dx. This particular choice of diffusion coefficient is exac % that required to 

damp the oscillations completely, while a smaller value would not do so. This coefficient 
also smears (diffuses) the solution badly, and reduces the overall accuracy of the approx- 
imation to first order. In the limit as Ax -+• 0 and At — ► 0, with a = aAt/Ax fixed, the 
dissipative term vanishes (as it must, for the finite difference approximation to remain 
consistent with the original inviscid equation), but it vanishes as a first order term. 

The quest for a good nonlinear scheme may therefore be considered as a quest for a 
good diffusion coefficient. We can write the equation to be solved as 


du du _ d ( du\ 
dt + a dx dx \ dx) ’ 


(2.38) 


with the understanding that k is nonzero only because our grid intervals axe nonzero. 
In the limit as Ax and At vanish, k must also vanish. Therefore we explicitly define 
K ex. At, in order that k vanish in the appropriate limit. (Why do we not define k cx Ax 
instead? Because to do so leads to confusion and the violation of an isotropy condition in 
multidimensional problems, as we’ll see later. In addition, while At always has units of 
time, the spatial grid parameter may have different units in different problems-angular 
units, for example, in problems with cylindrical geometry. Thus it would be difficult to 
come up with a general expression for k, which is proportional to the grid spacing.) We 
also require that the relative sizes of the convective (a du/dx) and dissipative terms be 
independent of the time step, as indeed they are in Eq. (2.37). These requirements are 

met for 2 

«c = kAt (2-39) 

M 

where k is a dimensionless, non-negative constant. Note that k is never less than zero, 
and is independent of the sign of a. The case k = 1/2 recovers the monotone first order 
result of (2.37). The case 0 < k < 1/2 reveals that while a linear scheme may not achieve 
the best of both worlds, it can achieve the worst: an oscillatory first order scheme! 

Note that the above scheme works just as well on a nonlinear equation, such as 
Burger’s equation, where a(u) = u. Then a, /c, and a = a(u)At/ Ax are functions of 
position, although At is not. In this case a du/dx should be written as df jdx, f = u / 2, 
in order to ensure conservation and the correct shock jumps (as in Chapter 1). 

Eq. (2.39) may be converted into a nonlinear scheme by defining 

k. = kAtj—M, (2.40) 

M 

where u depends on the local numerical solution U. We can define u in such a way as to 
satisfy u < 1 always, to have v « 1 near discontinuities and oscillatory regions, and to 
vanish as a first order quantity in regions where U is a smooth, continuous function. In 
finite difference form we write 


Vx (*^) i = - u i) ~ Ui- 1)] , (2.41) 
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for which 


*•'+1/2 — 2^ Ki + *»'+l)> 


Ki = kAt-^-Vi, 


Oi — Qj 


At 

Ax 


Of the many possibilities we could pick for i one of the most useful 


is 


Ui = 


I Ui+i ~ 2 Ut + U^ | 
Ui+i - Ui I + I Ui - Ui-! 


(2.42) 

(2.43) 

(2.44) 


This definition has = 1 if U { is a local maximum or minimum. Thus K i+ x / 2 will be 
a maximum only if U{ and f7,-+i are both local extremes, which is by definition a spurious 
oscillation to be damped. Note that the local wave speed a,- and Courant number o t - 
appear, so that the above prescription may be applied to the general case of spatially 
varying wave velocities. The quantity a?/ |a,-| oc |a,|; hence the diffusion coefficient is 
larger when the wave velocity is larger, which is as it should be, since the solution (and 
spurious oscillations) will evolve more quickly in regions where a,- is large. 

The constant k sets an upper limit to the diffusion coefficient, and is left as a free 
parameter to be set by the user. Different problems will give the best results for different 
k values. The choice k = 1/2 should eliminate spurious oscillations; however, in many 
cases k = 1/2 will be excessive, and a smaller value should be chosen. 

The preceding analysis has focused on the choice of a second order central difference 
approximation for the du/dx (or df/dx ) term. One could choose the fourth order 
approximation (2.11), in which case the above scheme would still reduce to a first order 
solution for u = 1, but which is not guaranteed to be monotonic. Nevertheless, one finds 
in practice that those oscillations which do occur are more easily damped (i.e., require 
smaller k values) than is the case when second order approximations are made. 


2.5 The Method of Lines in Two Dimensions 


The extension of the preceding techniques to two or more dimensions is straightforward. 
We begin with a single conservation law in two spatial dimensions: 


or 


du df dg 

ai + ai + a^ = 0 ’ / = / M' » = »(«). 


du t du du df , dg 

ha— + 6— = 0, a = — , b = — 
oy du du 


dt ~ dx 

The two equations are equivalent. If a and b are constant, the solution 


is 


(2.45) 

(2.46) 


, ^ / ax + by \ 

u(z,y,() = Uo -- t J, (2.47) 

which is a wave moving at constant speed, with velocity components a in the x direction 
and b in the y direction. 
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The second order semi-discrete approximation to Eq. (2.45) is 

dUjj _ fi+ij - Ihli g ' j+1 — g 'i~ l (2.48) 

dt ~ 2Ax 2Ay 

in the absence of dissipation. The grid points are assumed to be uniformly spaced, with 
x,- = tAx, y,- = jAy, and Ufi = C/(xj, y,-,t n ). 

Dissipation is introduced by adding diffusive terms to Eq. (2.45) in the form 


dt 


+ 


du , £/,££_ iL ( K §l\ + — ( K — \ 

dx + dy dx \ dx) + dy \ dy) ’ 


(2.49) 


where k — + 0 as At — » 0, as before. We define k. in terms of the local Courant number 
and the largest local velocity component c, = max(|a,|, |&,|)‘ 


(2.50) 


Kij = kAt^-Vij, ( 2 - 51 ) 

Oij 

where a useful (but not unique) definition for i is 

( \Uj+i i - 2 Ujj + l^ij+1 ~ 2 ^t i + Ujj - il \ (2.52) 

Vii ~ maX ( I C/,+1,- - Uii I + I Uii - Ui-xiY \U ij+ i - C /,,-1 + |c/, - c/,-1-1 1 ; ' 

The same definition for k is used in both dissipative terms in Eq. (2.49). This unique- 
ness in the definition of /c is required by the transformation properties of scalar fields, a 
subject which will be taken up in detail in Chapter 4. 

The second order semi-discrete approximation to Eq. (2.49) is then 


dUi y 

dt 


fi+lj — fi - lj _ 9*1+1 9tj-l 
~ 2Ax 2Ay 

+ [k.+1/2j(C/.+1j - Ua) - K i-l/2j{Uij - C/,-1,)] 

+— ^ [ K «j+ 1 / 2 ( t7 o+i “ U 'i) ~ ~ C/,-1)] > 

/C.-+1/2, = h*, + K ,+i/2 = 2 + 


(2.53) 


where n 

(2.54) 

X — 

The stability properties of this scheme are the same as for the one dimensional case of 
section 2.3, and with max <7, < <T mai , for o as defined by Eq. (2.50), and o max as given 

in section 2.3. 
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2.6 Exercises 

1. Using the fourth order approximation to df jdx\i given in Eq. (2.9), compute the 
sixth order approximation. 

2. Compute second and third order one sided approximations to df/dx | f , starting with 
the first order one sided approximation of Eq. (2.6). 

The following exercises are numerical solutions of the two initial value problems 

du df 

di + di = °' / = att ’ a=1 ’ 

u(x,0) = Ui(x) or u 2 (x), 

e -[(z-0.25)/0.l] J 

f 

1 x < 0.25, 

0 x > 0.25. 

In all cases let Ax = 0.01, x< = t'Ax, and compute the numerical solution U { for grid 
points i = 0 through 100. Define initial values for D, for t = -2 through 102 from the 
initial conditions, but do not recompute U_ 2 , U. u U W1 , or U W2 at the new time steps, 
(in other words, hold the boundary values fixed with time.) Compute 100 time steps 
at a Courant number o = 1/2. At time step 100, plot the analytic solution u x {x - at) 
or u 2 (x at) as a continuous unmarked curve, and the numerical solution U? n = 100 

as dots at the positions (x„C7"), « = 0, . . . , 100. Use the computer to plot Is well as’ 
compute the results. 

3. One sided (upwind) scheme (2.15) with u x and u 2 as initial data. 

4. Lax-Wendroff method (2.22) with u x and u 2 as initial data. 

5. Method of lines: scheme (2.31) with ( PU ), = -(/ i+1 - / 4 _ x )/ 2 Ax, using u, end Uj 
as initial data. 

6. Same as 3, but adding the dissipative term of (2.41) to the right hand side (i.e., to 

(PU)i as defined in problem 3), with k given by Eqs. (2.42)-(2.44), and k =0 0 2 
and 0.5. ’ ' ’ 

7. Same as 6, but using the fourth order approximation (2.11) to df/dx. 

8. Same as 3, but using Burger’s equation in conservative form (1.22) and the initial 
value problem of Exercise 3, Chapter 1. Compute the time step size at the beginning 
of each time step according to At = oAx/U max , where U max is the largest absolute 
value of the numerical solution (max/°° |U?|) on the grid at the beginning of the 
current time step. The left boundary condition must now be changed to one of 
antisymmetry, given by U-i = -U x , U_ 2 = —U 2 . Does the numerical solution 
match the analytical solution? 


ui(x) = 
u 2 (x) = 


Chapter 3 

Stability Analysis 


Chapter 2 presented several approximate finite difference solutions for one dimensional 
conservation laws. The criteria used to select an appropriate technique are stability, 
accuracy, and efficiency. Of these stability is the most important, as the formal accu- 
racy and efficiency of a method are irrelevant if the method is unstable in practice. The 
stability limits for the methods in Chapter 2 were stated without proof. In this chap- 
ter I will present a general technique for determining the stability of finite difference 
approximations for wave equations. 

3.1 The Consistency Condition and the Lax Equiv- 
alence Theorem 

The partial differential equations we wish to solve may be written in the general form 

^ = Pu, (3.1) 

dt 

where P is an operator which acts on u to give its time derivative. The quantity Pu may 
be nonlinear and contain any combination of powers of the spatial coordinates, time, 
or unknown elements u fc , or any spatial derivatives of these combinations, but may not 
contain time derivatives of u. 

As in Chapter 2, we consider only the pure initial value problem, for which the initial 
data are given, and the time dependent solution is to be obtained, over an infinite spatial 
domain with no boundaries or boundary conditions. 

Replacing the spatial derivatives in P with suitable finite difference approximations 
yields the semi-discrete approximation 

^ = (PU),, (3.2) 

at 

where U, is the approximate solution for U at the tth grid point. 

Given data at time level t n , we integrate Eq. (3.2) to time level f n+1 = t n + At. Thus 
the new value of U" +1 is some function of the old values Uj* +Jb , k = — oo,...,oo. The 
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functional relationship between the new and old values may be written 


U," +1 = C(At)U", 


(3.3) 


where the operator C(At) acts on U f - and depends on A t. 

We know from section (2.1) that the quantity (U," +1 - 
to the time derivative of U, , hence V 


U ,^) /At i 


is an approximation 


C(At)U" -U" 
At 

must be an approximation to PU. 

The consistency condition therefore requires that 


|C(At) - I 


At 


p u .(0 


0 as At — » 0, 


(3.4) 


where ||/|| is any valid norm of / and I is the identity operator [5j. 

The consistency condition looks imposing, but it is simply a formal statement of an 
intuitively meaningful concept namely that the finite difference approximations which 
are used m the numerical solution of a differential equation must yield that equation in 

the limit as all grid spacings go to zero (i.e., asAx^Oand A1-+0 for the methods of 
Chapter 2). 

The consistency condition is automatically satisfied by any one-to-one replacement of 
derivatives by valid finite difference approximations, so the insistence on consistency may 
seem redundant. However, we saw in Chapter 2 that dissipative terms were required in 
t e finite difference solution to discontinuous problems, even though these terms do not 
represent finite difference approximations to derivative terms appearing in the original 
differential equations. Thus an approximation containing these added dissipative terms 
must be defined m such a way that the dissipation vanishes in the limit of zero grid 
spacing in order to retain consistency (otherwise our finite difference method is solving a 
problem different from the one whose solution we want). The dissipative terms in section 
.4 vanish as At — ► 0, because k oc At . If we had defined k so that /c ^ 0 as At 0, 
Ax — ► 0, then the approximation would have been inconsistent, and guaranteed not to 
converge to the correct answer in this limit. 

The definition of consistency given above leads to the Lax Equivalence Theorem [5]: 

Given a properly posed initial-value problem and a finite difference approximation to it 

that satisfies the consistency condition, stability is the necessary and sufficient condition 
for convergence . 

The proof of the theorem assumes P to be a linear operator, but in practice the same 
result is observed to hold for nonlinear systems as well, provided the dissipation is chosen 
properly. 

Having established that stable and consistent schemes will converge to the correct 
answer, we proceed next to define and analyze the stability of finite difference approxi- 
mations. 
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3.2 The Von Neumann Method for Stability Analy- 
sis 

According to Eq. (3.3), the operator C(At) advances the solution from time t to time 
t + At. Thus initial data U? yield data at time t n = nAt through the repeated application 
ofC(At), 

U? = [C(Af)] n U?. (3.5) 

A scheme is stable if there exists a r > 0 such that the set of operators [C(At)] n is 
bounded for 0 < A t < t. In other words, the numerical solution may not grow without 
bound (provided that the correct solution does not). Conversely, a numerical scheme 
which is unstable will exhibit unbounded growth, and this growth is virtually always 
exponential. 

Note that computing a bound for [C(At)] n is not trivial, as C(At) is an operator, not 
a number. Fortunately, a straightforward technique for transforming the operator C'(Ai) 
into a number does exist, when C (At) is a linear operator. The technique is known as 
the Von Neumann method and is described below. 

4 . • i 

3.2.1 Fourier Analysis and the Linear Wave Equation 

We consider once again the linear wave equation 

— -|- a— = 0, a = constant, (3.6) 

dt dx 

and examine the behavior of the numerical solution for u. Fourier analysis turns out 
to be a convenient tool for this examination. The Fourier transform of a function / (x) 
is the frequency spectrum of the function, /(w), and is a function of the frequency lo. 
If / = f(x,t ) and we take the Fourier transform of the x dependence, the spectrum is 
f(oj, t ). If in turn f(u, t) is a monotonically increasing function of t for some frequency w 0 > 
f(x, t) will be an unbounded function, since the frequency component wo grows without 
bound. 

The Fourier transform of the x dependence of a function u(x,t) may be written 

/{u(x,t)> = -i= f e~ iwx u{x,t)dx = u{u,t), (3.7) 

v27T 00 

and is finite if u(x, t) -*■ 0 sufficiently quickly as x -» ±oo, which we will assume through- 
out. 

The derivative dujdx transforms according to 
3r {£ u(x ’ ,) } = 

which is obtained upon integration by parts. 


(3.8) 
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Taking the Fourier transform of the linear wave equation (3.6) gives 

du 

— + luau = 0. (3.9) 

Notice that Eq. (3.9) is an ordinary differential equation for any particular choice of 
w. The effect of the Fourier transform has been to replace a single partial differential 
equation with an infinite number of ordinary differential equations. 

The solution of Eq. (3.9) is 

u(u,t) = e~ iuat u((j, 0), ( 3 . 10 ) 

and therefore 

u(u>, t) 

u(w,0) (3-11) 

Although the phases of the components of the frequency spectrum change with time, 
their magnitudes do not. Indeed, this is the behavior we would expect, since we know 
the solution to be a traveling wave of constant shape. 

3.2.2 Fourier Analysis for Finite Difference Approximations 

Previously we have taken U ? to be the numerical approximation to the exact solution at 
the point (x,,t n ). For the remainder of this chapter, we will adopt a slightly different, 
but equivalent, interpretation, and assume that U* is a function defined by 

U i+k = U(x + kAx,t) , (3.12) 

which we happen to sample at discrete points. Thus any derivative approximation, such 
as 

— " ~ El±l Z Ebl - + - U{x - A x,t") 

dx i ~ 2Ax 2 Ax ’ (3.13) 

x i 

is also a continuous function, and we may compute its Fourier transform. 

Assuming Ax to be a fixed constant, we have 

7{U[x,t)} = U, (3. 14 ) 

TU(x + lcAx,t)} = e ,h *U, £ = u>Ax. (3.15) 

The variable £ is referred to as the dimensionless frequency, and is an angular quantity 
(hence dimensionless). 

j can now compute the Fourier transforms of finite difference approximations. Let 
*i 6 z m) Ui be an mth order approximation to dU/dx\i. Then we represent the one sided 
approximations by 


«W-t7, = Ui - U, tW*U, = U M - Ui, 


(3.16) 
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and the centered approximations by 

tf’l'i = - Oi-O. ( 317 ) 

S^U, = i|8 {U M - Vi-,) - (0i« - Oi-i)]. (3-18) 

The second derivative, d 2 U jdx 2 \i is approximated by ^ Ui , to order m, and we 

have 

Si = U i+1 - 2 Ui + Ui- 1 . (3.19) 

The transforms are 

= (1 - «“*)£, 7{«? )+ 0i} = (e i€ - l)£ , (3-20) 

J{<5(2)Z7,} = (t sin £)U , (3.21) 

7{6^Ui) = *^[8 sin £ — sin(2f)]17 , (3.22) 

O 

T{6 2 x W Ui} = -4sm 2 {Z/2)U. (3.23) 


3.2.3 Stability Condition for Numerical Methods 


We can now state the Von Neumann stability condition for numerical approximations to 
the linear wave equation. All one step linear finite difference solutions to the linear wave 
equation (3.6), implicit or explicit, can be written as 

oo oo 

E = E '*!£*, (3.24) 

k—~oo k= — oo 

where the l k and r k are constants. Taking the Fourier transform of Eq. (3.24) gives 





(3.25) 


from which we see 


U n+1 = g(t)u n , 



r k e 




i k e^r 


(3.26) 

(3.27) 


where 47(f) is known as the Fourier amplification factor. 

The Von Neumann stability condition is the following [5]: Stability of finite difference 
approximation (8.24) requires that 9(f), os defined by Eq. (3.2 7), satisfy 


| ff (f) | < 1, for all f in - tt < f < n. (3.28) 

Otherwise the numerical solution grows exponentially with the number of time steps. 1 

The preceding formalism may appear complicated, but is simple to implement in 
practice, as the following examples show. 

1 A more general condition allows for exponential growth when such growth is a valid solution, and 
requires g(f) < 1 + O(At) [5], This situation occurs when a source term bu is added to the right hand side 
of (3.6). For most problems condition (3.28) is sufficiently general, and will be used in these notes. 
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3.3 Simple Upwind Schemes 

Let us first consider the upwind schemes of section 2.2.1. The Fourier transform of 
Eq. (2.15) is 

U n+l = [1 - a{ 1 - (3.29) 


ff(0 = 1 — o(l — e **), (3.30) 

= 1 - a(l - cos £) - t'asin £. (3.31) 

Multiplying Eq. (3.31) by its complex conjugate gives the square of the absolute magni- 
tude of g as 

lff(0l* = 1 “ 4cr (l - <?) sin 2 (£/2). (3.32) 

The worst case has £ = n, for which 


Iff 001* = l-4a(l -a) < 1, (3.33) 

a(l - a) > 0, (3.34) 

and we must have 0 < o < 1 for approximation (2.15) to be stable. A similar calculation 
shows that the stability bounds for approximation (2.16) is -1 < o < 0. Thus the 
stability of the explicit upwind schemes depends on the sign of the velocity a. 

Consider next the implicit upwind scheme of Eq. (2.17). Its Fourier transform is 


U n+l = U n -a{l-€~ i( )U n+1 , 

ff = 1 - < 7(1 - e~ i( )g. 


9(0 = 


1 

1 + <7(1 — cos 0 + ursin 0 


(3.35) 

(3.36) 

(3.37) 


1 + 4a(l + cr) sin 2 (^/2) ’ 

Again the worst case occurs at £ = n, and the condition |ff(£)| < 1 implies 


(3.38) 


^(l+^)> 0, (3.39) 

and the method is stable for any a > 0. Similarly, approximation (2.18) is stable for any 

a < 0 . 


3.4 Simple Centered Methods 

The simplest centered method we can have is the approximation 
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The Fourier transform of (3.40) gives 

g(£) = 1 - iff sin £, ( 3 -4l) 

|<7(£)| 2 = 1 + ff 2 sin 2 £, (3.42) 

so that approximation (3.40) is unconditionally unstable! 

In Chapter 2 we found a dissipative term which, when added to the simple centered 
difference approximation, produced an upwind scheme which we now know to be stable 
for suitable time steps. The dissipative term is given in Eq. (2.37), and a first order time 
integration of (2.37) is 

= U? - i<T (t£, - «r_i) + Wh - W? + or.,). (3.43) 

The Fourier transform of (3.43) gives 

g(£) = 1 — |ff|(l — cos 0 - *ff sin £, (3-44) 

|ff(£)l* = 1 + 4 i° 2 - M) sin*(£/2), (3.45) 

and the stability condition is 

ff 2 — |ff| < 0, (3.46) 

which is satisfied for |<r | < 1, and is independent of the sign of the velocity, unlike the 
original one-sided schemes. 

We conclude that adding dissipation to a centered scheme which is originally unstable 
may produce a stable scheme. Thus we now consider a more general scheme of the form 

ur l = y, n - y>(DT + 1 - vi i) + i«(or + . - w + or.,). < 3 - 47 ) 

where a = nAt/Ax 2 and k is a constant diffusion coefficient, as in Eq. (2.38). The 
Fourier transform of Eq. (3.47) gives 

g = 1 — ia sin £ — 2cc(l — cos £), (3.48) 

| g (£)| 2 = 1 - 8asin 2 (£/2) + 16a 2 sin 4 (£/2) + ff 2 sin 2 £. (3.49) 

The condition |jr(£)| < 1 yields the equation 

/(£) = [2a - ff 2 + (ff 2 - 4a 2 ) sin 2 (£/2)] sin 2 (£/2) > 0 for all £. (3.50) 

Eq. (3.50) specifies a relationship between a and ff. Analyzing (3.50) in detail serves 
no purpose for the present, but for a given o it puts lower and upper limits on the range 
of a values which will lead to a stable approximation, and it defines a maximum value 
of |ff| above which any approximation will be unstable for any value of a. 
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3.5 The Lax-Wendroff and MacCormack Methods 

The Lax-Wendroff and MacCormack methods yield Eq. (2.25) when applied to the linear 
wave equation. The corresponding Fourier transform gives 

9{0 = 1 - *<7sin £ + < 7 2 (1 -cosf). (3.51) 

Using Eq. (3.50) with a = a 2 /2 gives 

(a 2 - a 4 ) sin 2 (£/2) > 0, (3.52) 

or |cr| < 1 for stability. These methods have the same stability limit as Eq. (3.43), 
although Eq. (3.43) is a first order accurate method, while the Lax-Wendroff and Mac- 
Cormack schemes are second order accurate. 


3.6 The Method of Lines 


In the previous sections we have analyzed the stability of a method by writing the new 
value of U , f/ t - , explicitly in terms of the old values U”. Obtaining such an expression 
for a method of lines integration such as (2.30) or (2.31) is a tedious business, due to 
the many substitutions (substeps) involved. It is simpler to take the Fourier transform 
of the semi-discrete form first, and then impose a time integration algorithm [9]. 

We’ll begin by considering the semi-discrete form of the linear wave equation, with 
added dissipation given by formulas (2.39) and (2.41), 

dU, a f \ k 

It = "A x 6 * u< + - 2U< + U ‘-J’ ( 3 - 53 ) 


k = kAt 



(3.54) 


where 6 ^ is an mth order undivided difference operator, as given in equations (3.16)- 
(3.18). Eq. (3.53) may be written 


dUi 1 
dt ~ At 

and its Fourier transform is 


1-oS^Ui + k\a\{U i+l - 2 Ut + U , ;_!)], 


dU 1 . 

dt ~ ~At kU ' 


A = — io-p( m )(£) — 4fc|a|sin 2 (^/2), 
while the functions are 


(3.55) 


(3.56) 

(3.57) 


P (2) (0 = sin(0, P (4) (0 = ^[8sin(£) - sin(20]- 


(3.58) 
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Re(A) 

Figure 3 . 1 : Stability regions for Runge-Kutta methods of order p. Copyright © 1971 by 
Academic Press, Inc. Reprinted by permission of Academic Press, Inc. 

At this point we have two options. The first is to perform the time integration of 
(3.56) with either method (2.30) or (2.31), both of which give 

9 — 1 + A- + 2^ 2 + + 24 (3.59) 

The condition |ff(£)| < 1 defines a region in the complex plane known as the stability 
region , such that the method is stable if A lies within the stability region. 

Much of the effort of obtaining a stability limit has so far been avoided by performing 
the Fourier transform before the time integration. Tlowever, determining the stability 
limit for A remains a non-trivial task, and a desire to avoid the effort involved leads us 
to the second option, which is to look up the stability limit of our time integration for 
Eq. (3.56) in a book on the numerical solution of ordinary differential equations. 

Page 120 of the book by Lapidus and Seinfeld [ 10 ], reproduced here as Figure 3.1, 
gives stability regions for several methods. The figure shows the stability regions in 
the complex plane of A for Runge-Kutta methods of order p, p = 1 , ... ,5. The p = 4 
curve applies for methods (2.30) and (2.31). This curve has imaginary intercepts at 
z = ±2 \/2i, and real intercepts at z = 0, —2.785. The symmetry of the curve about the 
real axis implies that the sign of the velocity a has no effect on the stability limits. 

If k = 0 , no dissipation is added and A = -ioP( m \£) is purely imaginary, and must 
lie in < A < \2\f2. Hence 


<xP (m) ( 0 l < 2\/2> 


(3.60) 
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and the maximum allowed value of a is determined by the largest value of -P^(f). 
We quickly determine that Pj£} x = 1, and Pf£} x — 1.3722, hence a < 2>/2 for stability 
if a second order approximation to du/d x is made, and a < 2.0612 if a fourth order 
approximation is made. (In both cases k = 0. Dissipation will decrease the stability 
limit even though it improves the quality of the solution.) 

At the other extreme we consider the purely dissipative problem, with a = a = 0, 
and replace k\a\ in (3.57) by tcAt/Ax 2 . Then A = — 4/c^^ sin z (^/2), and —2.785 < 
A < 0, from which we see that /cAf/Ax 2 < 0.6963 for stability. In real diffusion or heat 
conduction problems, where k is some nonzero function of the solution u, the maximum 
stable time step is proportional to the square of the grid spacing. Reducing the grid 
spacing by a factor of two decreases the allowed time step by a factor of four. A finely 
spaced mesh may therefore require a prohibitively small time step for an explicit scheme 
to be stable, which suggests that an implicit scheme might be more efficient, though 
considerably more difficult to implement. 

For the added dissipation we are considering here, /cAt/Ax 2 = k\o\, which vanishes 
as <7 — » 0, so we cannot obtain a purely dissipative problem while a is nonzero. However, 
if a is sufficiently small, we can make k large enough that the second term of Eq. (3.57) 
dominates over the first, and we may conclude that k is restricted by the requirement 
k\a\ < 0.6963 for stability. 

A general expression relating the maximum value of k as a function of the Courant 
number o probably cannot be derived analytically, but is not really necessary. If |cr| does 
not exceed roughly half of its upper limit, one can in practice make k large enough to 
damp any oscillations. The practical condition is to compute the time step from 

At = Ax-pj, (3.61) 

M 

where a is a Courant number which is less than the maximum value, and to set k to an 
optimum value obtained by experimentation (e.g., |a| = 1, A: = 0.3). 

3.7 Stability Analysis in Two Dimensions 

Stability analysis for two or more dimensions is very similar to the one dimensional case, 
and we can use the one dimensional results with only slight modifications. We start 
with the Fourier transform in two dimensions (x,y) of a function u(x,y,£), which may 
be defined 

1 roc roo 

7{u(x,y,t)} = — dx dy e~ x ^ >xX+WyV ^ u{x,y ,t) = u(to x ,uj v ,t). (3.62) 

Z71 J — oo J -oo 

It follows that 

7U[x + mAx,y + nAy,t )} = (3.63) 

where 

( = w,Ai, rj = u> y Ay. (3.64) 
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As in Eq. (3.24), any one step linear finite difference solution to the linear equation 
Eq. (2.46) may be written 


£ s u U&] it ,= £ r u V^ tJ+l . 

Jk,I=— oo k,l= — oo 


(3.65) 


where the s^i and r « are constants. 

We take the Fourier transform of Eq. (3.65) and define the amplification factor g(£, ij) 
by 

U n+1 =9{Z,v)U n ' (3-66) 


The Von Neumann condition is then 

|ff(6,*?)|<l, -7T < £ < 7T, -7T<»7<7r. 


(3.67) 


3.8 The Method of Lines in Two Dimensions 


As in the one dimensional case of section 3.6, we write the differential equation in semi- 
discrete from first, 


dUg 

dt 


a S lm) U xMj/.. 

Ai^ U,} A y ’ U 13 

—(«+« - 2K, + Vi-ij) + A(^y+i - 2C/„ + 0#-,), (3.68) 



c = max(|a|,|6|), 


with «j to be defined below, and then take its Fourier transform to get 


dU 

dt 



(3.69) 


(3.70) 


A = _4,c ^ sin2 ^/ 2 ) ~ 4,c ^ sinJ (^/ 2 )- ( 3 - 71 ) 

Once again we simply look up the stability region for A in the complex plane for the 
particular time integration method we use. If either method (2.30) or (2.31) is used, we 
get as given by Eq. (3.59), but with A as defined above in Eq. (3.71). 

Consider now the case of no dissipation ( k = 0). Then A is purely imaginary, and 
achieves its extreme values for r) = ±£ such that P^(£) is maximized. Since the extreme 
values of A on the imaginary axis for which the method is stable are ±t’2\/2, we see that 



A tPjZl < 2y/2. 


(3.72) 
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(The absolute value signs are used in (3.72) to cover the case where a and b have opposite 
signs, in which case the most extreme A occurs for P( m )(£) and P^(r]) maximum in 
magnitude but of opposite sign.) If we define the Courant number in two dimensions by 


a = 


'A + I^L 

Ax Ay 


) At, 


(3.73) 


we recover the familiar stability limit of (3.60), written as 

oPj£l < 2y/2. (3.74) 

If m = 2, 0 < 2v/2; if m = 4, <r < 2.0612. We select a stable value for a and compute At 

o 


from 


At = 


j£l + J*L‘ 

Az ' Ay 


The case where dissipation dominates leads to the result 


+ 


1 


Ax 2 A y 2 


/cAf < 0.6963, 


(3.75) 


(3.76) 


which for k as defined by Eq. (3.69) puts an upper limit on k for any given At, or vice 
versa: 




c 2 At 2 , 

k < 0.6963. 

a 


(3.77) 


3.9 Stability Limits for Fluid Dynamics Problems 


As we will see in Chapter 5, the three equations of one dimensional fluid dynamics may 
be transformed into an equivalent set of coupled nonlinear wave equations. The wave 
velocities are the characteristic velocities of the system, namely v — c, v, and v + c, 
where v is the fluid velocity, and c is the speed of sound. It is reasonable to expect that 
the largest characteristic velocity determines the time step for the whole system. We 
therefore pick a value for o < o maz and compute At from 


Ax 

A t = a t-t, (3.78) 

maxfCjJ 

c* = M+c, (3.79) 

where the maximum is over all points on the grid. 

In two dimensions the x direction characteristic velocities are v x — c, v x , v x + c, while 
the y direction velocities are v v — c, v v , v y + c. Hence we set 

At = (3.80) 


max 


( _L _£jl 
\Ax ^ Ay 


)' 


c x = \v x \ + C, C v =|v v | + c, (3.81) 

where the maximum is over all points on the two dimensional grid, and a < cr max . 

The stability properties of the dissipative terms are essentially the same as in the 
previous section. A more detailed discussion of dissipation for nonlinear systems will be 
given in Chapter 4. 
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3.10 Exercises 


1. Verify Eq. (3.11) in the one and two dimensional cases, using the explicit solutions 
for u(x,t) and u(x,y,t) given in Chapter 2. 

2. Derive Eqs. (3.20)-(3.22). 


3. Compute the stability limit for the approximation 


ur' =v?~ jw - 


Could you have predicted this result? 


4. Derive Eq. (3.50). 

5. Find the real and imaginary intercepts of the equation |gr (A) | = 1, with g given by 
Eq. (3.59), and A taken to be a general complex number. 





Chapter 4 

Artificial Viscosity and 
Conservation Laws 


The exercises in Chapter 2 demonstrated the need for artificial viscosity in finite difference 
approximations which represent discontinuous solutions. The artificial viscosity took the 
form of added dissipative, or diffusive, terms, which acted to damp numerical oscillations 
and spread out the discontinuity over a few grid intervals. 

One requirement which the added dissipation must satisfy is the consistency condi- 
tion. The dissipative terms must vanish in the limit as all grid intervals go to zero. Con- 
sistency is a necessary condition for valid artificial dissipation, but is far from sufficient. 
One could define an infinite number of terms which satisfy the consistency condition but 
which do not have the appropriate behavior. 

The purpose of this chapter is to define the Appropriate behavior” for artificial 
viscosity. However, we must first understand the general properties of the original non- 
dissipative system, and it is with the system of inviscid fluid dynamics equations that we 
begin. 


4.1 Conservation Laws and Tensor Fields 


Tensor calculus is the natural language in which to express conservation laws for field 
quantities because it allows one to write general equations which are coordinate invariant, 
i.e., which are valid in all coordinate systems. The simpler techniques of vector calculus 
suffice to describe the conservation of scalar fields, but are inadequate to describe the 
properties of vector fields. As the momentum density is a conserved vector field, the 
tensor description will be employed. A brief review of tensor calculus is given in the 
following sections, but a comprehensive derivation of the formulas is not possible in this 
text. The interested reader is referred to the book by Weinberg [ll], or any text on 
tensor calculus, for the details. 
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4.1.1 Tensors 


A tensor is represented symbolically by a symbol (such as a Greek or Roman letter) 
followed by some number of upper and lower indices (superscripts and subscripts), as 
in T ab C d e . The indices may be either letters or numbers. A numerical index refers to 
a particular tensor component, while a letter index refers to all components, with the 
understanding that giving the index a numerical value singles out a particular component. 

The total number of indices is the rank of the tensor. (The above example has a rank 
of 5.) In a coordinate system with n dimensions, each index represents a number from 
1 to n 1 , and a tensor of rank r has n r components. For example, the 9 components of 
the second rank tensor A £ in 3 dimensions are A}, A], A 3 , A \ , A\ , Af, A 3 , A 3 , and A®. 
Tensors of rank zero have no indices, and are called scalars. Tensors of rank one have 
one index, and are usually called vectors. Tensors of higher rank have no special names, 
and are referred to generically as tensors. Tensors of rank two, however, are conveniently 
written in matrix format, and are sometimes treated as matrices. 

Tensors are defined mathematically as objects which satisfy certain coordinate trans- 
formation properties. Suppose we have two coordinate systems, one of which has co- 
ordinates written as x a = (x 1 ,! 2 ,! 3 ), and the other of which has coordinates x a ' = 
(x 1 ,x 2 ,x 3 ). One example of two such systems are the rectangular coordinate system, 
in which x° = (x, y,x), and the spherical coordinate system, in which x°' = (r, 6,<f>). We 
can define the following coordinate transformation matrices to relate the two systems: 


A“,= 


dx a 

dx a> 



(4.1) 


A tensor Tg in the unprimed coordinate system may be transformed into an equivalent 
tensor in the primed system according to 

' T$ = K'<>T b % (4.2) 

where repeated indices, with one up and one down, are summed over their range (i.e., 1 
to n), so that Eq.(4.2) is a compact representation of 


Tg 


E E 

a=lfc=l 


(4.3) 


A tensor with upper indices only is contravariant. An example is the coordinate 
velocity vector u a = dx a /dt, which transforms according to 


u 



dx °' dx * _ A a'<*X« 

dx a dt 0 dt 


~ A“ u a . 


(4.4) 


A tensor with only lower indices is covariant. The gradient s a of a scalar S is covariant, 
since 


= dS_ _ dx^_dS_ dS 
Sa dx a> dx a> dx a a> dx a 


A“.s„. 


(4.5) 


1 Or 0 to n — 1 in general relativity. 


4.1. CONSERVATION LAWS AND TENSOR FIELDS 


43 


The coordinate transformations defined above assume that the origins of the original 
and transformed systems are the same. Translations (transformations which move the 
origin) are excluded. 2 The allowed transformations are (l) rotations of a particular 
coordinate geometry about some axis (e.g., a 45° rotation of a rectangular system about 
the axis x = y = z); (2) transformations from one coordinate geometry to another (e.g., 
rectangular to spherical); or (3) any combination of (l) and (2). 

For example, rectangular and spherical coordinate systems are related by 


x = rsinfl cos(/>, 

y — r sin 9 sin <f>, 
z = r cos 9, 


r = y/x 2 + y 2 + z 2 , 

tan9 = ^HZ, 

tan <f> = y/x. 


(4.6) 


4.1.2 Coordinate Systems and Metric Tensors 


Let t represent the time coordinate, and x° = (x 1 ,® 2 ,® 3 * * * ) represent the spatial coordi- 
nates. The vector x a is a contravariant vector. The coordinate system denoted by x° 
need not be rectangular coordinates; for example, the choice x a = (r, <f>,z) represents 
cylindrical coordinates, while x a = (r, 9,<f>) represents spherical coordinates. (Note that 
all indices run from 1 to n .) 

Now consider an object which is moving through our system. At any one time it 
has spatial coordinates x a , but the coordinates change with time. Consequently we may 
define a coordinate velocity u a by 

■* - (4 ' 7 > 


i.e., u 1 = dx l /dt y u 2 = dx 2 jdt , and u 3 = dx s /dt . The coordinate velocity is a contravari- 
ant vector whose components are the time derivatives of the spatial coordinates. The 
coordinate velocity is not the same as the physical velocity , represented throughout by 
v a , which is the rate of change of distance along the coordinate axes and which is not 
a tensor! For example, the coordinate velocity in spherical coordinates is while 

the physical velocity components are v r = f, v$ = r$ y v^ = rs\n6<j>, where the dot denotes 
derivatives with respect to time. 

Clearly we need to relate the coordinate velocity to the more intuitive notion of 
physical velocity. To do so, we must have some way to relate coordinate changes to 
the distances spanned by them. In tensor calculus, these two quantities are related 
by the metric tensor The metric tensor is a symmetric second rank tensor, whose 
components are usually functions of the coordinates. 

Suppose our test particle moves from position x a = (x 1 ,# 2 ,^ 3 ) to x a + dx a = (x 1 -f 
dx x ,x 2 + dx 2 ,x 3 + dx 3 ), where dx a is a set of small coordinate differentials. In doing so, 


2 The four dimensional spacetime of general relativity involves more general transformations, as time 

is one of the coordinates to be transformed. Lorentz transformations assume coincident origins for the 

original and transformed spacetime coordinate systems, which means that the spatial origins of the two 

systems coincide at time t = 0. Lorentz boosts are particular transformations which relate two systems 

moving with respect to each other, and thus produce time dependent translations. 


44 


CHAPTER 4. ARTIFICIAL VISCOSITY AND CONSERVATION LAWS 


the particle traverses a distance ds given by 

ds 2 = g a f,dx a dx b . 


(4.8) 


The symmetry of the metric tensor implies that only n(n + l)/2 of its n 2 compo- 
nents may be unique (six components in three dimension j). A general metric in three 
dimensions may be written schematically by a 3 x 3 array of numbers, similar to a matrix, 


( 



0n 

012 

013 

Qab — 

021 

022 

023 


^ 031 

032 

033 


Qab ~ 

9ba- 



(4.9) 


where by symmetry 

(4.10) 

In orthogonal coordinate systems, where the coordinate axes are perpendicular to 
each other, the metric tensor is diagonal, and may be written 


( 


Qab 


h\ 


h\ 


j 


(4.11) 


where hi, h^, and h$ are the scale factors for the coordinate directions (and the superscript 
2 represents the square of the number, not a tensor index). 

A diagonal metric simplifies tensor calculations. For example, we need to know the 
contravariant (inverse) metric g ai , defined by 




(4.12) 


where S b is the Kronecker delta, 


*1 


1 a = b, 
0 a b. 

The inverse of a diagonal metric g a b is simply 


( 


9* b = 


hi 


-2 


hi 2 


K 2 


(4.13) 


(4.14) 


The inverse of a general tensor is considerably more complicated. 
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The metric tensor is of crucial importance, not only because it defines the distance 
relationship for a coordinate system, but because it is also used to transform a contravari- 
ant index into a covariant index, and vice-versa, a process referred to as “raising and 
lowering indices.” If a is a lower index on a tensor, and we wish to raise it, we multiply 
by g ab to get a new tensor of the same name with an upper b index where the old lower 
a index was. Similarly, multiplication of a tensor with an upper a index by g a b produces 
a new tensor of the same name with a lower b index where the old upper a index was. 
Some examples are 

u b =g ab u a , u h = g ab u b , A ab cd = g ae g b fA efed , (4.15) 


and so forth. 

We can now compute the physical velocity, defined above as the rate of change of 
distance, per time, in the coordinate directions of an orthogonal coordinate system (de- 
scribed by the diagonal metric of (4.11)). Let the physical velocity have components 3 
(t 7 i, t> 2 , V 3 ). The distance traversed by a motion which changes the jth coordinate by dx J 
and does not affect the others is 

ds = hjdx J (not summed), (4-16) 


hence 

Vj = hjU* — Uj/hj (not summed). (4-17) 

Thus in spherical coordinates the contravariant (coordinate) velocity u a is (r,9,<f>), the 
covariant velocity u a is (r,r 2 0 ,r 2 sin 2 04 >), and the physical velocity v a is (f,r 0 ,r sin #<£). 

The metrics and inverse metrics for the three most often used coordinate systems are 
given below. 


Rectangular x a = (x, y, z) 



1 



\ 

1 

Sab — 

1 


_a& 

> 9 = 

1 







Cylindrical x a = (r, 4 > , z) 




\ 


\ 

1 

Sab — 

r ! 


ab 

. 9 = 

r " 2 







(4.18) 


(4.19) 


3 Remember, the physical velocity is not a tensor! 
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Spherical x a = (r, 




\ 




\ 


l 




1 

9ab — 

r 2 



Q ab = 

r- 2 




r 2 sin 2 0 j 




r -2 sin -2 0 j 


4.1.3 Combining Tensors 

New tensors may be generated from old tensors in a variety of ways. Two tensors of the 
same type may be added to give a third of that type: 

K = K + Cl. (4.21) 

The product of two tensors of rank m and n is a new tensor of rank m + n, 

Aab cd = B ai c cd , ( 4 . 22 ) 

provided no two indices are the same. If two indices (one up, one down) are the same, 
the implied summation produces a new tensor whose rank is two less than the original, 
as in 

Rl = S ac bc , (4.23) 

T\ e = F a de G bc d \ (4.24) 

The summation process is often called contraction, and is a frequent occurrence in 
tensor calculus. As the number of indices involved increases, the number of terms in 
the sum increases rapidly. For example, the contraction of a tensor R abcd with itself to 
produce a scalar R, according to 

R = RabcdR aicd , (4.25) 

has 3 4 = 81 terms in 3 dimensions. 4 


4.1.4 Covariant Differentiation 


Another way to generate new tensors is to differentiate old ones. The derivative of a 
tensor with respect to a scalar is a tensor of the same type as the first, as in Eq. (4.7) 
for the coordinate velocity. 5 

The derivative of a tensor with respect to a coordinate direction is usually not a 
tensor. Such an object is represented either by the usual derivative notation, or by a 
subscript comma followed by an index indicating the direction of differentiation, as in 


*/ , ^l_ r c 

dx° La ’ dx c ~ ab 


( 4 . 26 ) 


4 Such an operation occurs in 4 dimensions in the calculation of tidal forces near black holes, and the 
corresponding sum has 256 terms! 

6 In relativistic calculations, the time coordinate is part of the four dimensional coordinate vector, not 
a scalar, and a different definition for velocity is used. 
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Note that the gradient ( djdx a ) of a scalar is always a tensor (a covariant vector), while 
gradients of tensors of rank one (vectors) or higher are in general not tensors. 

One can define a more general derivative operation than the above which does yield 
tensor results when applied to a tensor (i.e., which is coordinate invariant). The operator 
is known as the covariant derivative and it reduces to the ordinary derivative of (4.26) 
when the original tensor is a scalar, or when the coordinate system is rectangular. 

The definition of the covariant derivative involves the connection coefficients T£ c , 


•po a dp 

1 be ~ 9 1 dbc, 

(4.27) 

r 1 ( dg d b dgdc dg bc \ 

dbe 2 y dx c dx k dx d ) 

(4.28) 

Note that is not itself a tensor, although for computational convenience I have used 

the metric to raise and lower its first index. The connection coefficients for rectangular, 
cylindrical, and spherical coordinate systems are given in Appendix A. It is a convenient 
fact that the connection coefficients for rectangular coordinates are all zero. Note that 
the coefficients are symmetric on the last two indices (r£ c = r® t ), due to the symmetry 
of the metric. 

The covariant derivative of a tensor in coordinate direction a is denoted by the sub- 
script ; o. The following are useful examples: 

e 

ii 

ft 

(4.29) 

B a . b = B a , b + Y a ci B\ 

(4.30) 

B a<b = Ba,b — KbBc, 

(4.31) 

c ab , c = c ab tC + v a dc c dl + r^c ad , 

(4.32) 

✓-t a pa i s~id pd 

v s b t c — ^ b t c i 1 bc^ 

(4.33) 

Cab.c = C'ofc.c - Y d c C db - T d c C ad . 

(4.34) 


Covariant derivatives exist for tensors of all ranks, but it will not be necessary to illustrate 
any more of the possibilities than given above. 

One useful property of the definition of covariant differentiation is that the covariant 
derivative of the metric is always zero: 

9ab;e = 9*; c = 0. (4-35) 

Thus the metric commutes with the covariant derivative operator, and we may write 
expressions such as 

T ab c - e = ( g cd T abd ). e = g cd T abd , e . (4.36) 

Covariant derivatives also obey the usual sum and product rules as ordinary derivatives. 
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The covariant derivative as defined above adds a single covariant (lower) index to a 
tensor. We can define a “contravariant derivative” by raising the derivative index with 
the metric: 

T ab c' e = g de T ai c , d = (g de T ab e ). d . (4.37) 

The contraction of a covariant derivative with an index of the original tensor is often 
called the covariant divergence. Examples of covariant divergence are V a ]a , T ab ,b , and 
so forth. While the preceding formulas apply just as well to the divergence as to single 
derivatives, there are some special relations which simplify the divergence calculation. 

If we write the metric g ab as a matrix, then we can compute its determinant. Let 

g = vl det(<7 a {,)|. (4.38) 

Then it can be shown that 

_ a 1 dg d 

r *“ = ?a? = a? Ins ’ < 4 - 39 > 

and the divergence of a vector V a is 


V.^ = -~(gV a ). 


;a 


g dx a 


(4.40) 


4.2 Fluids as Tensor Fields 

Three types of fields are encountered in inviscid fluid dynamics. The density (p), pressure 
(p), and total energy density (e) are scalar fields. The coordinate velocity (u°), momen- 
turq. density (m u = pu a = pg ab u b ), mass flux (pu°), and total energy flux ([« + Pi**) 
are vector (first rank tensor) fields. The momentum flux (pu a u b + p6 b ) is a second rank 
tensor. 6 

The scalar conservation laws for density and total energy are 

+ (/>“*);« = o, (4.41) 

+ [(« + p)«°];o - o. (4.42) 

The momentum equations are more complicated. First we define the (symmetrical) 
momentum flux tensor [12] 

T ab = pu a u b + pg ab , (4.43) 

or, equivalently, 

T b = pu a u b + pS b . (4.44) 

The momentum equations are then 

dm a m i 

-gf- + T a b ., b = 0, (4.45) 

®If it weren^ for the momentum flux, we could get by with the simpler formalism of vector calculus, 
but such, alas, is not the case. 
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dm a 

~df 


+ (m 0 u 6 + p5*);6 = 0. 


(4.46) 


The system is not complete without the following auxiliary relations. Let v be the 
velocity magnitude; then 

v 2 = u a u a = g ai u a u b . (4.47) 

For a perfect gas, the pressure is given by the equation of state 

P = h ~ 


‘ = e ~2 PV 


(4.48) 

(4.49) 


where e is the thermal energy density. 

The above equations are useful for formal analysis, but their numerical solution re- 
quires writing the equations explicitly in terms of ordinary derivatives. Using the diver- 
gence relations of the preceding sections gives the more familiar-looking results: 

(4.50) 

(4.51) 

(4.52) 


dp Id , 

at + ^ > = °’ 

l4a? Is(e + pK, = 0 ’ 


dm. 


+ — 


i a 


(gm a u h ) + 


dp 


dg 


cd 


dx • + 2 pu ‘ Ud ax‘ 


= 0. 


dt ' g dx b 

The last term in (4.52) is a “centrifugal force term,” which arises in non-rectangular 
geometries. 

Integral forms of the conservation laws may be obtained by multiplying each equation 
by the volume element gdx 1 dx 2 dx 3 and integrating over the volume desired. If the 
velocity goes to zero at the boundary of the volume, we get 

f f J P 9 dx Y dx 2 dx s = 0, (4.53) 

^ J J J e g dx x dx 2 dx 3 = 0, (4.54) 

and the total mass and energy of the system are constant. The momentum integral is 


7tl U m -<" lx ' dx ’ dx ’ = -I U{ 


dp 1 dg 

aT- + 2 ^Tx 


ed\ 


gdx l dx 2 dx 3 , (4.55) 


so that the individual components of momentum density are not conserved in general, 
even though the momentum field as a whole is. Note, however, that the rectangular 
components of momentum density are conserved, since in rectangular coordinates g = 1 
and dg cd /dx a = 0, so that 

J J J m a g dx l dx 2 dx 3 = 0 (rectangular coordinates) (4.56) 

assuming p is constant over the boundary. 

The reader who is impatient to see what the complete set of fluid equations looks like 
may turn to Appendix A, where the equations are written out in full for the three most 
popular coordinate systems. 
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4.3 Dissipation for Tensor Fields 

Let U be a generalized density (scalar, vector, or tensor) and F = F(U) be a generalized 
flux, satisfying the conservation equation 


dV 

dt 


+ div(F) = 0, 


(4.57) 


where div(F) is the divergence operator defined in the previous sections. Eq. (4.57) is 
representative of all the fluid equations. 

If the relationship F = F(U) is nonlinear, shock waves are likely to arise in the 
solution of Eq. (4.57). Even if the flux function is linear, we may still have contact 
discontinuities in the solution. 

Therefore the numerical solution to Eq. (4.57) is likely to require the addition of 
artificial viscosity, or dissipation, to prevent nonphysical oscillations and nonlinear insta- 
bilities in the presence of discontinuities. 

The added dissipation must be formulated with care, in such a way as to produce 
physically meaningful solutions. One could define many dissipative terms which would 
damp oscillations but introduce nonphysical behavior into the solution. 

We add dissipation to the numerical solution by replacing Eq. (4.57) with the equation 


au 

dt 


+ div(F) = D, 


(4.58) 


where D is the dissipative term. For the presence of D in Eq. (4.58) to give meaningful 
numerical solutions to the invtscid Eq. (4.57), we require the following conditions: 

L ® does not contribute as a source term for conserved fields when integrated over 
any finite volume, although it may contribute a flux across the boundary of that 
volume; 

2. D acts to diffuse (spread out) sharp gradients and damp oscillations; 

3. D is written in a general form which is valid in all coordinate systems; 

4» D is isotropic— there are no preferred coordinate directions; 

5.. D vanishes as all grid intervals (At, Ax°) go to zero (consistency condition). 

The first point requires that D be written as the divergence of some kind of flux. The 
second point requires that the flux be a diffusive one, involving derivatives of order no 
higher than first. The third point is satisfied automatically if D is written in tensor form. 
The fourth point requires that the diffusion coefficient which appears in the flux be a 
scalar, not a tensor of nonzero rank. The fifth point requires that D be proportional to 
some power of a grid interval, and for the fourth point to be satisfied that grid interval 
must be At, and not one of the spatial intervals Ai a , 
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In Chapter 2 we found that a monotonic (oscillation free) solution for single wave 
equations could be obtained by adding to the right side of the equation the term 

/ 1 <7 1 Ax 2 \ dul 

\ 2 At J dx \ ’ 

which is a diffusive term with diffusion coefficient However, we now require that 

Ax not appear explicitly in the coefficient, so as to preserve isotropy in the multidimen- 
sional case where we will have multiple grid intervals Ax°. We write D as a diffusive 
term of the form 

D = div[/cgrad(U)] (4.59) 

where k is a function of position. Eliminating Ax from the above diffusion coefficient in 
favor of cr, a, and At gives 

/c = kAt^-, (4.60) 

M 

where a is the Courant number. This definition will be extended to multidimensional 
systems in a later section. 


d_ 

dx 


4.3.1 Scalar Fields 

If U is a scalar field U, the flux F is a vector with components F a , and the conservation 
equation with dissipation may be written 



(4.61) 

(4.62) 

(4.63) 


4.3.2 Vector Fields 

If U is a covariant vector field 7 U„, and the flux F is a tensor F*, the conservation 
equation with dissipation is 

^T + F ‘*={ KU ?)y ( 4 - 64 ) 

Working out the covariant derivatives on the right and making use of relation (4.39) gives 


dt a ’ g dx c 


be 

99 * 


dUa -T d ab U d 


dx b 




(4.65) 


7 Similar results apply to contravariant vector fields, but the momentum density is most conveniently 
expressed in terms of covariant fields, so we consider only the covariant case here. 
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Unfortunately we cannot eliminate the connection coefficients r? c in favor of simpler 
expressions, as we did for the scalar equation. The presence of the connection coefficients 
produces the inhomogeneous term on the right, which implies that the components of 
U a are not conserved in general, except when rectangular coordinates are chosen. In 
rectangular geometry all r?,. = 0, and we see that the rectangular components of U a are 
conserved. The dissipative operator conserves precisely those components of the field 
as does the original momentum equation! That the same conservation properties hold 
in both cases is not an accident, but is a basic property of the geometry of coordinate 
systems. 


4.4 More on the Dissipation Coefficient 

The fluid equations with dissipation are given in full in Appendix A, and will not be 
repeated here. Now we will look at the definition of the dissipation coefficient ac in more 
detail. 

As shown in Chapter 2, the choice 


a 


k = kAtj-r, a = a(U) 


At 

Ax’ 


(4.66) 


reproduces the monotone first order scheme of (2.37) when k = 1/2, but does this 
definition (and its generalization to higher dimensions) have the conservation, isotropy 
and coordinate invariance properties we seek? 

The answer is a qualified yes. Isotropy can always be maintained by defining ac 
uniquely at each grid point, so that dissipation occurs at the same “rate” in all coordinate 
directions. Conservation is guaranteed by the divergence form of the dissipative term, 
so only coordinate invariance needs to be investigated further. 

The definition of ac given above is not coordinate invariant, nor should it be. The 
dissipation depends on the choice of grid, as it must in order to satisfy the consistency 
condition, and will change as the grid parameters change. Thus ac is not invariant under 
grid changes in any given geometry. Moreover, coordinate rotations in multidimensional 
problems will alter the characteristic velocity components, changing a and cr, and hence 
k. Thus ac is not invariant under coordinate rotations, and is not a scalar. 

What happens, then, when we do perform a coordinate transformation on some phys- 
ical problem and solve the problem numerically in both coordinate systems? Will the 
two solutions agree? They will, provided both have adequate resolution. By “agree” 
I do not mean agree exactly (an obvious impossibility, since the two solutions will not 
even be computed at the same grid points), but the solutions should agree to within 1% 
or so, except possibly at discontinuities. The diffusion coefficient ac adjusts automati- 
cally to changes in the grid and coordinate system so as to apply the correct amount of 
dissipation needed in each case. 

Although strict coordinate invariance is neither necessary nor desirable, strict isotropy 
is both. If ac is computed differently in the diffusive terms for different coordinate di- 
rections, then it is really “simulating” a tensor ac?. For any given problem the results 
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may appear reasonable. However, to reproduce those results in a^new coordinate system 
would require transforming ac£ into = A“ Ajr/Cj, and using kJ in the new system. If 
the transformation is not performed, the two solutions will in general not agree. If the 
transformation is performed, the new is likely to contain off diagonal elements and 
be an explicit function of the transformed coordinates, even though the original K a h was 
a diagonal tensor with constant elements. 

Tensor diffusion coefficients may have a role to play in the real world (such as diffusion 
and heat conduction in anisotropic materials), but they are of no use in artificial viscosity. 


4.5 Dissipation in One Dimension 

This section illustrates one method for raising the order of the dissipative terms to provide 
a globally second order accurate solution, for single equations and the one dimensional 
fluid equations. To get a second order scheme, simply multiply the definition of (4.66) 
by a function v[U) which is also first order, 


a 


k = kAt — i/, 
a 


(4.67) 


where 1 / w 1 in oscillatory regions and near discontinuities, and is a first order quantity 
otherwise. One definition which has proved quite useful is 




1 Uj+i - 2Uj + Ui - il 
\u i+1 - Ui\ + \Ui - Ui.,1' 


(4.68) 


where Ui is the numerical solution at grid point x,-. 

In summary, we make a finite difference approximation to the dissipative equation 


du 

dt 


+ 


df d ( du\ . 


(4.69) 


T 

(4.70) 


as in Chapter 2. The semi-discrete approximation is 

^ - Ui) - Ki-HliUi - U^)}, 

dt Ax Ax A 

using the definition (4.67) for /c, where the wave velocity a = df /du, and the constant 
k is constrained by fc|cr| < 0.6963, for the integration schemes (2.30) and (2.31). (The 
choice k = 0.3 is often effective.) Note that the quantities /, k, a, and a - aAt/Ax are 
local functions of U and will generally vary from place to place. 

Now consider the more difficult problem of fluid dynamics in one dimension, choosing 
spherical geometry for the coordinates. Retaining only the radial terms from Eqs. (A.43)- 
(A.45) gives 

dp . 1 d , 2 1 d f_ 2 dp\ 

pu ) 


(A *71} 
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where 


dm r 

dt 


dt Id,,, 

T t + 7*d-r [r < e + P > U l 

1 d . 2 dp 

+ 7 , T r (r m ' u > + * = 



(4.72) 

(4.73) 


u r = m r /p, ( 4 . 74 ) 

P = (7 - 1) (e - ^Pu r2 ) = (7 - 1) (e - ^m r u r ) , (4.75) 

and 7 is a constant (5/3 for a monatomic gas). Here u r is the coordinate velocity, which 
happens to be the same as the physical velocity v r in this case. 

The semi-discrete approximation to these equations is 


dpi 

dt 


de f - 

dt 


dm r i 

dt 




r?Ar 


r 2 Ar 2 [ r * +1 / 2 /c *+ 1 /2 (p*+i Pi) r i— 1/2 — l /2 (p« ~ A'-i)] > 


r?Ar 2 [ r * ? + 1 /2 AC *+ 1 / 2 ( e *+ 1 “ e ') “ r i-i/2 K i-i/2i e i ~ «i-i)] , 
7 ^ 7 < 5 r m) (r, 2 m r ,u, r ) - 4z S r m) Pi 


Ar 


(4.76) 


(4.77) 


r ?Ar 2 [ r ‘ 2 + 1 / 2 ' c ‘ +1 / 2 ( m r*'+l ~ m ri) - r i-i/ 2 K i-i/ 2 (m r i - m r ,-i)] (4.78) 


— j/c,m r i. 

r i 


To compute the diffusion coefficients we need to define an effective wave speed a. In 
the case of a single equation, a was just the local value of df j du. In the current problem 
we have three equations, and there are three characteristic velocities (as will be shown in 
Chapter 5). The three characteristic velocities are u r — c, u r , and u r -f c, where c = \fypj~p 
is the speed of sound. The amount of dissipation needed is governed by the most rapidly 
changing wave. Hence we define the quantity c r to be the magnitude of the greatest wave 
velocity in the r direction, 


and set 


c r « — K| + c,-, 

C r > 


Oi = At 


Ar’ 


— ^ri» 

Then the dissipation coefficient is given by 


(4.79) 

(4.80) 

(4.81) 


a. 


Ki = kAt-r-^rVi, 

w 


K i+l/2 — 2^ K * /C ‘+ 1 )’ 


(4.82) 
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where i/< may be made a function of any solution component ( p , m, e, or some combina- 
tion). The choice u = u{p) has been found to be effective, so we define 


_ jPt+l 2 P» d~ Pi-ll 
_ \Pi+l - Pi\ + I Pi - P<-1 


(4.83) 


Note that the equations are undefined at r = 0, because of the coordinate singularity. 
The time derivatives become infinite at r = 0 unless u r (r = 0,t) = 0, in which case 


L’Hopital’s rule gives 


dp_ 

dt 


d . r . 

+ 3 ^(pu) = 


dt 

dt 


+ 3 ^[(e + p)« r ] 



(4.84) 

(4.85) 


m r = 


0 


r=0 


(4.86) 


We know from the condition of spherical symmetry that p , e, and p are symmetric 
functions of r, while m r and u r are antisymmetric, and can use these symmetry conditions 
to evaluate the finite difference approximations at and near r - 0. 


4.6 Dissipation in Two and Three Dimensions 

The simplest two dimensional problem is the single conservation law of Eqs. (2.45)- 
to which we add dissipation as follows: 

du df dg _ d ( du\ d_ ( 

dt + dx + dy dx \ dx) dy \ dy ) 


(2.46), 

(4.87) 


The i direction characteristic velocity is a = df/du, and the y direction velocity is 
b = dg/du. 

The semi-discrete approximation to (4.68) is 


dUjj 

dt 


1 fH/ 

Ax ^ 1 f ' 3 Ay v 9,3 


f }{Ui+ 1 j - Uij) - Ki- 1/2 j{Uij Ui-1 ,)] 

+ ~ Uij ^ - K 'i- 1 / 2 ( Ui i - 


(4.88) 


The local Courant number is 



Now select the dissipation on the basis of the maximum velocity component, 

Cij = max ( | a,,- 1 , |6 t y |), ( 4 - 90 ) 
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where 


^iy — h At j 

c *+i/2 j — 2 ( K, ‘i K *+i i) » /c *j-hi/2 = ~ (^t'j + 


(4.91) 


f ^ ’ - 2D » ± Bzi -I |Pil« - 2 tig + Pi,.,! \ ' 

VW'+i i - d.,1 + IK,- - £1,-1 ,ntl, i+ . - u,j] + \U„ - £/,,_. I J • 

Finally we consider the fluid equations in three dimensions, and in cylindrical coor- 
dinates (r, <p,z). The equations, with dissipation, are 

m + + + 

= i± (r„?l\ U. ± « /'„«£>') , a / . . 


= i A , 1 A , 5 / dp\ 

r dr V dr ) + r 2 d(j> ( d<j>) + dz \ K dz ) ’ 

dt + ^f r(e + p)u ^ + ^[( e + p) u 1 + ^[( c + p)« 1 

= — fricliU 


rSrC"ilr) T ,.S^^ T S^j' 

5 m r 15 , 9 d fin 

~dt~ + 7^ (rmrUr) + a^( m - u ^ + ^( m ^) + ^-rpu^ 2 


- IA/A,A^ , 1 ^ [ l dm ' 5 / dm r \ 

- rar[ rK ar J + ‘ ~]J + 51 (* 17 j «.*) 

1 /'dm, \ 


+ + £(»*«') + JW) + 


' ‘ d 4 > v 
dm* 

dr r 


1 d r 

+ T7T U 


r 2 d(fr { ^ d<f> 

. d i 1 (dm r 5m*\ 

M A)’ 

+ ~(rrn.«') + ±(rn.u *)+£(".,»*) + g 

= IA A d m A , 1A f d ( dn 

r dr V 5r J + r 2 (* d<£ J + dz ( ‘ K d 


+ rm T 


+ al *-57 


u r = m r /p, u+ = m*/(r 2 p), u' = m t /p, 

P = (7-1) e- ip(u r2 + r 2 u^ 2 + u z2 )] 

= (7 - i) « - ^ (m r u r + m,u* + , 


(4.97) 


(4.98) 


(4.99) 


(4.100) 
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and t is a constant (5/3 for a monatomic gas). The coordinate velocity «• is used in these 
equations, rather than the physical velocity v a (= h a u a , not summed), as the equations 
have a simpler (and less singular) form when coordinate velocities are used. Note that m* 

is an angular momentum, not a linear momentum: m * = pu * = g^pu* = r 2 p4> = rpv *. 
The semi-discrete approximation to these equations at the grid points (riy<f>j>Zk), 

where r,- = iAr, <f>j — 3&4>, %k — k&.z, is then 


dpi, k 

dt 


dcjjk 

dt 


dm r 


ijk _ 


dt 


-y^rt m) {riPi jk u r ijk ) - ^ m) (p, } kut jk ) - -^ s t { , m) (Pijk u ijk) 

+ rZr * [ r<+1 / 2K<+1 / 2 > k ( pi+1 > k ~ pijk ) ~ r *'- 1 / 2 ' c *'- 1 / 2 i k ( pi ’ k ~ pi ~ l **)] 

+ [*0+1/2 k{pxi+l k ~ Pijk) ~ *0-1/2 k{Pijk ~ Pij-i *)] (4.101) 

+ ~Kz* [ K ’> fc+1 / 2 ^ ,,Jk+1 _ Pijk ^ ~ Ki i k ~ 1 / 2 ( pi ’ k ~ ’ 

-^~5 r (m) [r,(eo* + PokKjt] “ +PH k ) u i}k} 

~ si m \{*ijk + P ijk ) uU \ (4,102) 

-1 fr, + x/ 2 *«+l /2 jfc(e.+l jk — «0'fc) _ r t-l/ 2 *i-l /2 jk{ e ijk ~ c «-l jk)] 

r»A r 2 1 

~*" r 2 A </> 2 k ~ e ' jk ^ ~ K,J_1 / 2 — C * J -1 *)] 

+ — ^ [ K ^ fc+1 / 2 ^’ jk+1 ” c< >*) _ *»jfc-i/2( e Ofc "" c 0*-i)] » 

^—<5 r (m) (r,m r i jk u r ijk ) - -rr4 m) ( mr ~ Az^ m ^ mr 


r,Ar’ r v ,JR ' A<£ 


+ 7 - [r <+ i/2*,+l/2 J*( m r .+1 jk - m r ijfc ) - r,-i/2*,-l/2 jk(»™ r O’fc m ' <-1 

r^Ar 2 L 

• V 


+ 


p+i t - wtr o* _ ™-<t> 0+1 k + rn* ijk 


( H2l- 

*0+1/2 k ^ A( ^ 2 r< 


Qk^ 


r i A<^ _ 

/ m r jjjfc - m T jj.. i t _ rrij, iik + p-i t \1 
-*0-i/2 fc ^ A< ^ 2r, JJ 

+ — ^ [ /C ^' fc+1 /2( m '- O'fc+l “ m r 0*) _ *Ok-l/2( m r 0* “ m r O'k-l)] 

1 f m, 0+1 k - m* 0-1 * A .. \ 


( 4 - 
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ijk 

dt 




A<£ 

1 






r i'+l/2«»+l/2 jk 


r,Ar 

~ r «-l/2 /c i'-l/2 ;fc 


«'+l jk in# jjk m<f> »+l jt + TTIq ijit 

A r 


2 r i+l/2 

/ m * Hk ~ ,-i jk _ , yt + ,_ x jt V 

V Ar 2r,— 1/2 J 


*, y+i /2 * ( — <J ’ +1 * >yt + r,- r i'tl * ± m - <&' 

V &<P 2 


r?A^ 

+ A? [*«*+•/» ( m * «*+i ~ m * •/*) - Kijt-i/ifmt ijt - m t 
k~n 

2A <t> 2Ar 


(4.104) 


+— /c f y A f ! 7lr ji±l * ^ mr 0-1 * __ m <t> «+! jk - jk j 


dm, 

dt 


r.Ar ^ ^ r,m * } ( m * <;*«£-*) ~ — <5, (m) (m* ,y*u* yt ) 


^ rfm) 

_ z; 4i R '‘ 


(4.105) 


1 r 

+ r,Ar 2 L r,+1 / 2AC ’' +1 /2 J*( m * »+i y* - m« «*) - r.-i/2«.-i/2 ,*(m, ijk - m, yjk )J 
+ r? A<£ 2 t*'** 1 / 2 *( m * </+* k~rn z i)k ) - /c,y_ 1/2 *(m, <yJt - m t t y_! fc )j 

+ Az 2 [' c ’2 fc+1 / 2 ( m * <y*+i ~ m * »y*) - «.y*-i/2 (m z ,- iJt - m, iyi -i)] . 

This system has the characteristic velocity components u r - c, u r , u r + c in the r 
direction; ru* - c, r u * ruf + c in the <f> direction; and u z - c, u*, and u* + c in the z 
direction, where c = is the speed of sound. To compute the dissipation, we first 

define the largest characteristic velocity magnitude in each direction according to 


define the Courant number as 


a = At(-^ + 


n + c > c* = \u e \ + c; 

(4.106) 

C <t> - c z 

rA(f> + Az ’ 

(4.107) 


and select the largest of the three values in (4.106) as the effective velocity a: 


a = max(c r ,c^,c*). 


(4.108) 
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The dissipation coefficient is then 


a ijk 

Kijk = kAt — 
&ijk 


(4.109) 


(4.110) 


where Vijk is taken to be 

( \pi+i jk - 2 Pijk + Pi-i jk\ 1p«j+i * ~ 2 P «J* + Eiizl *1 1 

v ijk - max y— jk -p ijk \ + \Pijk-Pi-i jk\' \Pij+i k ~ Pijk\ + I Pijk ~ Pi: -i k\ 

1 P» J fc + 1 %Pijk ~t~ Pijk - ll \ 

|Pij'fc+l — Pijk I "h I Pijk Pijk-l \J 

where p is the pressure. As usual, the dissipation coefficients at the half points are given 
by 

= i(/C,j* + K,+l jk), Kij+ 1/2 fc = «( K »'jfc + K ij+ 1 *)’ K ')k+ 1/2 — O K,,fc+1 ^' 

2 (4.111) 


i+1/2 jk ^ 


4.T Coordinate Singularities 

The fluid equations given above for cylindrical coordinates contain a coordinate singu- 
larity at r = 0. In the one dimensional spherical case we dealt with the singularity at 
r = o by requiring that the divergence term in each equation be finite, which allowed us 
to use L’Hopital’s rule to get a non-singular limiting form at the origin. We generally 
cannot require that each divergence term in an equation be finite at a singularity when 
dealing with two or three dimensional problems, as it is only necessary that the sum 
of all divergence terms be finite. (This problem becomes particularly acute in spherical 
coordinates, which is singular at r = 0 for all 6 values, and at 6 = 0, & = x for all r 
values.) It is quite possible that any given term may tend to infinity as at the singularity, 
even though the sum of all such terms in the equation is finite. 

The completely general three dimensional problem in non-rectangular coordinates 
is probably best handled by reverting to the rectangular coordinate definition for the 
divergence at the singular points, if individual terms are known to become infinite. If 
on the other hand particular symmetries of the problem dictate that individual terms 
remain finite at singularities, then L’Hopital’s rule may be used to obtain a limiting form 
for the equations. 


4.8 Dissipation for General Hyperbolic Problems 

So far we have looked only at the fluid dynamics equations. However, the formalism 
used above may be stated in very general terms for any first order system of hyperbolic 
conservation laws, in any number of spatial dimensions. In this section I will describe 
how to construct dissipative terms for an arbitrary system in three spatial dimensions, 
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representing conservation laws for conserved scalar and covariant vector fields The 

reader who is interested in less (or more!) than three dimensions should have no trouble 
making the appropriate changes. 

Our system will have i conservation equations for the scalar densities am am 

0 ( 0 , mvolvmg 1 flux vectors / (!) « where the subscripts are names fm the’ 

different densities, not tensor indices. 8 The equations are of the form 

dct(i\ 

dt + /( 0 ;<* = °- ( 4 . 112 ) 

We also have m conservation laws for the covariant vector densities Qi 

involving m flux tensors F m J, F m .> Note that tensor indict run 

trom 1 to 3, so there are 3m such equations, of the form 

+ ^ * 0 

Qt T r (') a u * 


Using the definitions for covariant derivatives gives the scalar equations as 

da(,) Id/ \ 

dt + gdx a ( 9 ^^~®’ , = 

and the vector equations as 
d(3 li)a Id/ .\ 

dt g dx b Mo. ) ~ ^abF(i)c =0, a = 1, 2, 3, t = 1, . . . , m. 


(4.113) 

(4.114) 

(4.115) 


The total number of differential equations is n = / + 3m. (There may be, and usually 
are some number of algebraic equations required to close the system, such as the equation 

of state.) It is convenient to define a single unknown variable vector U and the three 
associated flux vectors F„ by 


U = 


( 


t , \ 

a (l) 


fw 

a (i) 


fm 

£(1)1 


jpa 

^(i)i 

£(1)2 

, F a = 

F (l)2 

£(1)3 


^(1)3 

£(m) 1 


1 

£(to)2 


E pa 

^(m)2 

^ £(m)3 j 


jpa 

s 7 


(4.116) 


8 For the fluid problems discussed so far, / = 2, 
a (i) = p and a (2 ) = e. 

9 The fluid problem has m = 1, with £ (1 ) a = m a 
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so that the complete set of conservation laws may be written 

dV dFi dF 2 , 3F S 


dt 


+ 


dx 1 dx 2 ^ dx 3 


+ 


+ D = 0 , 


(4.117) 


where D is a vector of inhomogeneous terms. (By vector in this instance I mean any set 
of variables in the sense of linear algebra, not a tensor. The objects U, F a , and D do 
not satisfy tensor coordinate transformations.) 

The equations may also be written 


^ + A |U + b|^ + C^ + D = 0. 

dt dx 1 dx 2 dx 3 


(4.118) 


where A, B, and C are the n x n Jacobian coefficient matrices with components 


. dF xi 

Aij dUj ’ 


J3 Wi 
~ dUj ’ 


r dF >< 
dUj 


(4.119) 


Now let the eigenvalues of A, B, and C be A,, m, and ru, i = 1 , . . . , n. These eigenval- 
ues are the characteristic velocity components for the associated coordinate directions. 
A, for the x 1 direction, for the x 2 direction, and 17 ,- for the x 3 direction. The extreme 

values are 

Ci = max |A,j, c 2 = max |/x,j, C 3 = max |rj,j, 
and the largest of these is 

a = max(ci,c 2 ,c 3 ). 

The dissipation coefficient is therefore 


(4.120) 

(4.121) 


a 

k = kAt — is, 
o 


where the Courant number is 




Cl +~^ + 


\x 3 ) ' 


,Ax l Ax 2 Ax 3 

The quantity u at grid point (xj,xj,x£) is then 


(4.122) 


(4.123) 


u ijk = max 


1 /i+l jk ~ 2 /tjffc + /t -1 jk\ l/t; + l k tfjjk + /<j-l fc 


I*' * * * J n : - — r , 

|/»+l jk — fijk 1 + I fijk — ft- 1 J*l l/ij+l k ~ f'i k I \f*} k ~ f'j~ l fc l 
~ fjjk- 1| \ 

h I fijk ~ fijk - 1 1 / 


\fijk+l ~ fijk \ + 


(4.124) 


where / is any scalar variable which reliably indicates the presence of discontinuities and 
oscillations, such as pressure in the fluid case. 

The equations with dissipation are then 


1 d / a \ 1 d ( 

-aT I + ^(^) = ia^r 9 a* j 


(4.125) 
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for the scalar fields, and 


dt 


+ 


1 d 
g dx b 
1 d 
g dx c 


(wW) - 

»*»*■ 


(4.126) 


Kc*g bc 




for the vector fields, as in Eqs. (4.63) and (4.65). The approximation to semi-discrete 
form is handled as in the previous sections. 


4.9 Conservation Properties of Finite Difference Equa- 
tions 


The fluid equations may be written in a variety of forms, such as the conservative forms 
described in this chapter, and the characteristic and primitive forms given in Chapter 
5* Finite difference approximations may be made to the equations in any form, with 
the same formal accuracy, yet all approximations in this text have been made to the 
conservative forms. The reader may wonder if there is an advantage to be had in using 
the conservative form in numerical work. The answer is yes. 

First of all, the artificial viscosity terms require the conservative equations for their 
definition. However , once the terms are defined, one could write the equations in an- 
other form, so the necessity for artificial viscosity does not require that the numerical 
approximations be in conservative form, though they are more simply written in that 
form. 

More important is the notion of numerical conservation laws . The numerical ap- 
proximation to the conservative equations possesses an exact conservation law which is 
analogous to the original conservation law embodied by the differential equations, and it 
is for this reason that the conservative form is preferred in numerical work. 

A field u is said to be conserved if it satisfies the equation 


du 1 d , 

dt + g dx i ) “ °’ 


(4.127) 


is 


where /* is the flux of u in the ith coordinate direction. The conservation property i 
seen most clearly by integrating Eq. (4.127) over the volume a 1 < x 1 < b\ a 2 < x 2 < b 2 
a 3 < x s < 6 3 : ’ 

\du Id 
dt ^ g dx' 


JU'*d*d* f + ~(s/ i ) 


= 0 , 


(4.128) 


-J fJu S dx'dz' d x> = -J" fjgf') 

n b l b 7 ,42 .41 

, (9f 2 ) dx'dx s - (gf) 

1 a 3 */ a 3 J a 1 


dx 2 dx s 


a‘ 

b 3 


dx l dx 2 . 


(4.129) 
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The rate of change of the volume integral of u is given by the integral of the fluxes over 
the surface bounding that volume. 

A properly chosen numerical approximation to Eq. (4.127) retains an analogous con- 
servation property [13]. If the semi-discrete approximation is written as 

dUjjk 1 ^t+1/2 j k ~ — 1/2 j k j I ^t'j+l/2 k ~~ 1/2 k 

dt g ijk Ax 1 gijk Ax 2 

_l_ 1 Jjjk+1[2 ~ - 0 7 (4.130) 

Qijk Ax 3 


where F l = gf\ then its summation over the volume elements g^k Ax x Ax 2 Ax s gives 

i if V iit 9ii k Ax'Ax’Ax 3 E (^rVi/2 #. “ *?/, 

dt t,j,k=i y.*=i 

- E (*&«/> * - *?./. t)Ax 1 Ax 3 - £ (4.131) 

»,fc=i ».i =1 

The interior fluxes cancel out of the sums, leaving only the sums of normal fluxes 
over the boundary surface. The numerical scheme of (4.130) yields an exact numerical 
conservation law. A numerical approximation to a nonconservative equation does not 
produce a conserved numerical solution. 

Eq. (4.130) makes use of finite difference approximations of the form 

5F Fj+ 1/2 ~ -ft-i /2 (4.132) 

dx Ax ’ 

which involve mid-point fluxes EJ+i/i, rather than grid point fluxes 2*). The dissipative 
terms naturally have mid-point fluxes, since the dissipative flux is F = gK du/ dx, for 
which jj _ jj 

=«*+./.*>■/. A* (4 ' 133) 

to second order in Ax. The convective terms have grid point fluxes F, = gif,, however, 
from which the mid-point fluxes F i+ 1 / 2 must be constructed for the proof above to hold. 
The choice 


Fi+ 1/2 = + Fi) 


(4.134) 


recovers the familiar second order result, while 


F.+i/l = + Fi) - ^(Fi+, + Fi-,) (4.135) 

recovers the fourth order result. Zalesak [14] gives flux expressions up through eighth 
order. 

Note that it is not necessary to compute i^+ 1/2 explicitly, as given above, to ensure 
numerical conservation. The finite difference formulas previously used are properly con- 
servative, as demonstrated by the existence of mid-point fluxes such as (4.133), (4.134), 
and (4.135), along with the numerical conservation law (4.131). 
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For example, consider the single equation 


du Id 

dt + gdx^ 9 ^ ~ 



(4.136) 


where k is the artificial dissipation coefficient. The standard semi-discrete approximation, 
using a fourth order formula for the convective flux, is 


dUi 

dt 


1 


;[8(ffi+l/i+l fft-l/i-l) — (<7»+2/«+2 ~ Qi-lfi-l)] 


+ 


g< Ax 2 

which is exactly equivalent to 


gi 12 Ax 1 

[*<+l/2(tf.+l - Ui) - Ki-i/ 2 {Ui - t/,-— i)J, 


(4.137) 


dUi _ ffi+1/2 ~ Fj- 1/2 

dt Ax 


(4.138) 


where 

Fi+i/i = — (gf.+i/.+i + gift) - —{gi+tfi+2 + gi-ifi-i) - «,+ 1/2 , (4.139) 

and is therefore numerically conservative. 


4.10 Test Problems 


We now have at our disposal all the tools needed to obtain the numerical solution of 
some typical fluid dynamics problems. For the sake of simplicity, we consider only one 
dimensional problems for which the boundary values may be held constant with time, 
postponing the discussion of nontrivial boundary conditions until Chapter 5. Both rect- 
angular and spherical geometries will be considered, with spatial coordinates x and r, 
respectively. 

These one dimensional test problems are evaluated on a grid of unit length, divided 
into 100 subintervals (7 = 100, Ax = Ar = 0.01). The value of 7 in the equation of state 
is 5/3 throughout. The figures show the analytical solutions (solid lines) and numerical 
solutions (dots) for the density p, the pressure p, the momentum density m, and the 
velocity v, at a time t. The Courant number a , 


a = 


Ax ’ 


(4.140) 


is set to 1 throughout, where c x = max(|v z | + c) over the grid (and similarly for c r ). 

The first example is the shock tube problem, frequently used as a test for hydrody- 
namical codes (as in Sod [24]), and whose solution is given by Thompson [25]. At time 
i = 0 the system consists of two spatially constant, stationary states, adjoining at x = 0. 
The left state (x < 0) has p = p = 1, while the right state has p = 0.125, p = 0.1. 
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As time progresses, a rarefaction wave forms and moves to the left, while the contact 
discontinuity and a shock wave move to the right. 

Figure 4.1 shows the solution at time step 60, t = 0.265, with k = 0.3, just before 
any waves reach the boundaries. 

The shock wave is located at x = 0.5 and is resolved in « 3 grid intervals. Small 
oscillations occur around the shock. These oscillations could be eliminated by increasing 
k, at the expense of spreading out the shock (and other features in the solution). The 
contact discontinuity at i = 0.23 is spread out over m 10 grid intervals. Unlike shock 
waves, which keep a constant width, contact discontinuities in the numerical solution 
widen steadily with time. The nonlinear nature of the shock wave produces a tendency 
to steepen which is countered by the artificial viscosity, so the shock width remains 
constant. However, the contact discontinuity is a discontinuity in a linear wave, which 
has no steepening tendencies, and therefore spreads monotonically with time as it is 
diffused by the artificial viscosity. 

A general technique for obtaining sharply resolved contact discontinuities has yet to 
be formulated, and it is beyond the scope of this text to describe any of the numerous 
experimental methods which have been developed. As a first attempt, one may attempt 
to detect the presence of such a discontinuity and decrease the artificial viscosity there 
to reduce the spreading rate. However, some viscosity is needed to prevent oscillations 
near the contact discontinuity, so the spreading cannot be elimin&ted by this technique, 
although it can be slowed. It should also be mentioned that such discontinuities will 
spread even in the absence of artificial viscosity, as well as producing oscillations. 

The next problem considered is the Sedov solution for a spherically symmetric explo- 
sion, described by Landau and Lifshitz [12]. (N.B. Eq. (99.10) of the reference should 
have i/ 6 = 2/(7 - 2).) An amount of energy E is deposited at the origin at time t = 0. 
The resultant explosion produces a self similar solution bounded by an outgoing spherical 
shock wave. The density and pressure curves are sharply peaked at the shock, falling 
off rapidly for decreasing r. The density goes to zero at the origin, while the pressure 
flattens out and becomes constant with r away from the shock. The velocity is linear 
near the origin, but steepens somewhat near the shock. The similarity solution is valid 
as long as the shock is very strong, so that the density jump across the shock achieves 
its maximum value. 

The numerical solution is produced by distributing an amount of energy E = 1 over 
the innermost five grid intervals, in the form of thermal energy. Let the sum of the 
volume elements for the first five points be 

s 

vol = 4nArJ2 r l, (4.141) 

i=0 

then the initial pressure is 

Pi = (7 -l)E/vol, (4.142) 

for i 0, . . . , 5. For t > 5, p,' = 10 ®, We also have p — 1 and u = 0 everywhere. 
The resultant numerical solution is shown in Figure 4.2 at step 365, t = 0 614 and with 
k = 0.35. 
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The shock wave is very strong, and is resolved in two grid intervals. Once again 
a small amount of ringing is observed near the shock, which could be eliminated by 
increasing k at the expense of spreading out the shock further. The post-shock density 
should peak exactly at a value of 4 for this problem; the numerical post-shock density 
is » 3.25, which is low by » 20%. The density profile behind the shock is so sharply 
peaked that the numerical solution cannot change quickly enough to achieve the correct 
density peak. Increasing the resolution of the calculation by increasing the number of 
grid points would reduce the error in the peak height. However, the correct shock velocity 
is achieved even though the exact post-shock solution is not attained. 

The exact solution to the spherical blast wave has the density going to zero at the 
origin, while the pressure remains finite and the temperature goes to infinity. The numer- 
ical solution cannot reproduce this singular behavior at the origin. Instead the numerical 
solution has a very small but nonzero density in the center, and a high but finite tem- 
perature. The pressure near the center is accurately represented. Note that the sound 
speed is therefore very large at the center, and is in fact the largest velocity present in 
the calculation. The allowed time step is therefore determined by the speed of sound at 
the center, and not by the post-shock velocity, which is much smaller. 
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4.11 Exercises 

1. Let x° = (x,y,z) be rectangular coordinates, and let x a = ( r,9,<f> ) be spherical 
coordinates. Using the relationships of Eqs. (4.6), compute the transformations A“ ( 
and A*'. 

2. Using the results of problem 1, write down the coordinate velocities u° = (x,y,z) 
in terms of u a> = (f,#,<^), and vice versa. Do the results match what you obtain 
by taking the time derivatives of Eqs. (4.6)? 

3. Let x a = (x, y,z) be rectangular coordinates, and let x a = ( x',y',z ') be a different 
rectangular coordinate system. The z and z' axes are the same, while the angle from 
the x axis to the x' axis is 9. Compute the transformations A“, and A“ . Use these 
transformations to express the velocity x° = (x,y, z) in terms of x a> — (x',y',z'), 
and vice versa. Are the results what you expect? 

4. Verify that the metric g a b and its inverse g ai satisfy relationship (4.12), using metrics 
(4.11) and (4.14). 

5. Compute the connection coefficients for cylindrical coordinates. 

6. Write T a b ]e in terms of T a b and its ordinary derivatives, in cylindrical coordinates. 

7. Verify Eq. (4.35) for diagonal metrics. 

8. Derive Eq. (4.40) for V“. 

9. If the spherical components of the momentum density are not conserved, what 
meaning does “conservation of momentum” have? 

10. Show that the choice of (4.60) for the diffusion coefficient leads to the monotonic 
first order scheme of Eq. (2.37). 

11. Show that the mid-point fluxes (4.134) and (4.135) reproduce the standard second 
and fourth order finite difference approximations when used in the formula 

d-F _ Fj+ 1/2 ~ F^-i /2 
dx Ax 



Chapter 5 

Nonreflecting Boundary Conditions 


Numerical solutions to hyperbolic systems of differential equations, such as the fluid 
dynamics equations, are usually obtained over a finite region. The previous chapters 
have focused on techniques for obtaining the solution interior to the boundaries of the 
system. However, the time evolution of the system is governed not only by the state in 
the interior of the region, but also by the choice of boundary conditions. 

Many types of boundary conditions are possible. Fluid flows past solid objects (such 
as airplanes) require u solid wall” boundary conditions at the fluid-object interface. The 
normal velocity component must be zero at such a boundary, while the tangential velocity 
may or may not be zero (it is zero if viscosity is important in the boundary layer). Solid 
wall boundary conditions are discussed elsewhere (e.g., Beam and Warming [15], and 
Shih [16]), and will not be covered here. 

One common boundary type is the free boundary, which corresponds to no solid wall 
or other physical interface, and across which matter and information are free to pass. 
In this case the evolution of the system depends on waves which enter the system from 
outside the boundary. The outgoing waves are described by characteristic equations, 
while the incoming waves depend on information which is external to the model volume. 
It often happens that the external solution is not known, in which case one may choose 
a boundary condition which allows waves in the interior (including shock waves) to pass 
out of the interior without generating reflections. These nonreflectmg boundary condi- 
tions , originally postulated by Iledstrom [17] for the one dimensional case in rectangular 
geometry, have been generalized to the multidimensional case in arbitrary coordinate 
systems by Thompson [18], and are discussed below. 


5,1 Waves in One Dimension 


Consider first the one dimensional case in orthogonal (but not necessarily rectangular) 
coordinates. We have a system of n equations describing the behavior of n dependent 
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variables. Let U be the vector 1 of conservative variables, satisfying 


aU dF . 
~dt + ai + c= 0. 


(5.1) 


w here F is the flux vector, and C 1 is an inhomogeneous term not containing derivatives 
(which often arises from divergence terms in nonrectangular geometry). Eq. (5.1) de- 
scribes the conservation properties of the system; that is, it relates the rate of change 
of the integral of a field over a small volume to the flux of that field across the volume 
boundaries. 

An alternate form for Eq. (5.1) is the primitive system, with a vector of dependent 
variables U, which satisfies 

au au „ 

dt + A dx + C “ °’ ( 5 - 2 ) 

where A is an n x n matrix. The choice of primitive variable vector U is not unique 
(although the choice of conserved variables is), and could be defined as the conservative 
vector. The following analysis assumes that U and U are distinct. 

The two systems are related by 


with 


au au 

dt ~ p at ’ 

(5.3) 

aF _ au 

dx dx ’ 

(5.4) 

p dUi 

' 3 dUj ’ 

(5.5) 

Qy = 3F ‘ 

W l dUj ’ 

(5.6) 

<y 

PH 

1 

PH 

II 

<1 

(5.7) 

o 

II 

1 

0 

(5.8) 


and 


where P and Q are also n x n matrices. 

Now let 1,- and r, be the set of left (row) and right (column) eigenvectors of A, 
satisfying 

1.A = A,l„ ( 5 . 9 ) 

Ar,- = A,r„ (5.10) 

where the A, are n eigenvalues of A, ordered so that A! < A 2 < • • • < A„. (The system is 
hyperbolic if the eigenvalues of A are real.) Then we obtain a diagonal matrix A by the 
similarity transformation 


SAS -1 = A, 


(5.11) 


J In this chapter, the term “vector” refers to any set of variables, and is not a vector in the tensor sense 
described in Chapter 4. 
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where the rows of S are the left eigenvectors 1,-, the columns of S are the right eigen- 
vectors i\, and the matrix A is diagonal, with A* = A,-. (Note that the transformation 
follows from the orthogonality of the normalized left and right eigenvectors: ljr ; - = £«>•) 
Multiplying Eq. (5.2) by S gives 


s ™ +AS f + SC = 0, 
dt dx 


(5.12) 

‘ f 


or 




(5.13) 


in component form. Eq. (5.12) is the characteristic equation corresponding to the original 

forms (5.1) and (5.2). 

If we can define a new function V by 


dVi = 1,-dU + 1 iCdt, 


(5.14) 


then (5.13) becomes 


dVi 

dt 


+ A, 


dV± 

dx 


= 0 , 


(5.15) 


which is a set of wave equations for waves with characteristic velocities A,-, as in Chapter 
1. Each wave amplitude is constant along the curve C,* in the xt plane defined by 
dx/dt = A,. 

However, the definition of (5.14) generally can be made only if A and C are con- 
stant everywhere, or if no more than two differentials appear on the right side of (5.14). 
Otherwise the coefficients in (5.14) must satisfy Pfaff’s condition for the integrability of 
differential forms for the functions V% to exist [l], a condition not met for the fluid equa- 
tions. Nevertheless, the characteristic form of (5.13) holds true independent of (5.14). 


5.2 Nonreflecting Boundary Conditions in One Di- 
mension 

In the one dimensional case we wish to solve Eq. (5.1) over the region a < x < b. 
The problem is an initial boundary value problem, because both initial data in the 
region a < x < b and time dependent boundary conditions at x = a, b are needed for 
the problem to be well posed. Difficulty arises in the boundary condition specification 
because Eq. (5.1) generally contains eigenvalues of both signs at the boundaries, implying 
that waves are propagating into and out of the domain. It is therefore more fruitful to 
work with the characteristic form at the boundaries, since we can then consider each 
wave separately. 

The outgoing waves (those with A,- < 0 at x = a, and A, > 0 at x — b ) depend only 
on information at and within the boundaries. Thus those equations in the form of (5.13) 
which represent outgoing waves can be solved as is, or in any equivalent form. Properly 
designed numerical approximations to (5.13) for outgoing waves, which depend on one- 
sided finite difference approximations involving only interior and boundary points, will 
therefore be stable. 
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The incoming waves (with A, > 0 at x = a, and A* < 0 at x = 6) are another matter. 
They depend on data exterior to the boundary, and numerical approximations to (5.13) 
not involving exterior data will be unstable. Thus we need to know something about 
the exterior solution in order to specify useful boundary conditions. In some problems, 
particularly steady state aerodynamics problems, the far field solutions are known to a 
good approximation, and the appropriate values can be specified [19], [20]. 

For time dependent problems (and many steady state problems as well) it is often 
desirable to use so-called nonreflecting or radiation boundary conditions, which have the 
property of minimizing reflections from outgoing waves. A number of authors have at- 
tempted to create boundary conditions which have this nonreflecting property. Bayliss 
and Turkel [21] formulated a perturbation approach in which the perturbations about 
the desired steady state were expressed in terms of waves. They then imposed boundary 
conditions which annihilated the outgoing waves (i.e., prevented the generation of incom- 
ing waves). Engquist and Majda [22], [23] developed nonlocal, nonreflecting boundary 
conditions for linear systems. From their nonlocal conditions they derived a sequence of 
partially absorbing local conditions. Hedstrom [17] developed a nonreflecting boundary 
condition for the one dimensional rectangular, nonlinear case. As the only nonlinear 
condition, Hedstrom’s is by far the most useful for time dependent problems. It will be 
generalized to multidimensional problems and non-rectangular coordinate systems below. 

Hedstrom’s nonreflecting boundary condition [17] can be stated in the following way: 
the amplitudes of the incoming waves are constant, in time, at the boundaries. This is 
the same as saying that there are no incoming waves, as it is the change in amplitude 
which indicates a wave. Mathematically, this condition is 


dVi 


dt 


z=a t b 


= o, 


in terms of the wave amplitude V, of (5.14), or 


(5.16) 




= 0, 

x=a t b 


(5.17) 


in general, for those waves whose characteristic velocities are directed inward at the 
boundary. Eqs. (5.17) for the incoming waves and (5.13) for the outgoing waves com- 
pletely determine the solution at the boundaries. 

Note that Eq. (5.17) will not give the desired behavior for any problem which should in 
fact contain incoming waves. In such a case one must be able to specify something about 
the incoming waves. Fortunately, Eq. (5.17) seems to be adequate for many problems of 
interest. 

It is easy to write a general equation which automatically reduces to Eq. (5.17) or 
Eq. (5.13) for incoming and outgoing waves. The general form is 


f. dV „ , 

li^T" + + I 

i at 


* c ) 


z—a } b 


= 0 , 


(5.18) 
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where 




x.i.flS 

A » I « dx 


0 


for outgoing waves, 
for incoming waves. 


(5.19) 


Thus the characteristic and nonreflecting boundary conditions can be combined in a very 
natural way, unlike other extrapolation methods. . - 

The set of equations (5.18) is solved by a method of lines approach, in a way similar 
to (and along with) that for the conservative equations in the interior (as described in 
section 6). The function values are obtained at the discrete coordinate positions £<, where 


%i = a + * Ax, 


(5.20) 


Ai = (b-a)/I. ( 5 - 21 ) 

Eq. (5.1) is solved at the interior points, defined by 0 < t < I, while the boundary 
equation (in the form of Eq. (5.25) below) is solved at the boundary points, defined by 
i < o or t > I. (An interior scheme which uses fourth order finite difference operators 
requires two boundary points at each boundary.) The spatial derivatives in (5.19) are 
evaluated using one sided difference approximations 


au 

dx i 


J-(U, + 1-U 0 .'<0, 

Ax' 

— (u, — u,-x) »>/. 


(5.22) 

(5.23) 


To get an equation for the conservative variables, we first define £ as the column 
vector whose components are , and write 


S^ + £ + SC = 0, 
ot 


(5.24) 


which leads to 


^ + P(S- 1 £ + C) =0. 
ot 


(5.25) 


5.3 Waves in Two Dimensions 

In two dimensions the conservative system is 


au aF 

dt ^ dx ^ 


a g 

dy 


+ C'+C[ 


= 0 , 


(5.26) 


with C!,. and C' y representing non-derivative terms, as before. (Only the sum of the C 
terms matters; the sum has been partitioned into two terms to retain consistency with 
the one dimensional case.) We have the relations 


au au 

dt ~ r dt ’ 


(5.27) 
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dF au 

dG dU 

dx ~ Q dx ’ 

dy dy 

a = p- 1 q, 

B = P -1 R, 

Cx=P - 1 C' X , 

Cy = P- 1 ^, 


which relate the conservative form of (5.26) to the following primitive form: 


dV A dV dU - 

+ A— + B— + C x + C„ = 0. 


at 


dx 


dy 


(5.28) 

(5.29) 

(5.30) 


(5.31) 


Now let 1,-, r,-, and A,- be the left and right eigenvectors and eigenvalues of A. Sim- 
ilarly, let mi, s,-, and ^ be the left and right eigenvectors and eigenvalues of B. Then 
the matrices A and B can be put in the diagonal forms A and 1VI by the similarity 
transformat ions 

SAS" 1 = A, TBT -1 = M. (5.32) 

The rows of S (T) are the left eigenvectors 1, (m,), the columns of S _1 (T -1 ) are the 
right eigenvectors r,- (sj), and A (M) is the diagonal matrix of eigenvalues A,- (^). Then 
Eq. (5.31) can be rewritten as 


d U 

at 


+ S _1 AS 


au 

dx 


+ r‘MT^ 

dy 


+ C x + Cy = 0 . 


(5.33) 


which is as close to the characteristic form in one dimension as we can come unless S and 
T are the same (i.e., unless A and B are simultaneously diagonalizable), and which will 
be referred to as a characteristic form due to the presence of the diagonal characteristic 
velocity matrices. 


5.4 Nonreflecting Boundary Conditions in Two Di- 
mensions 

The two dimensional problem allows for an arbitrary number of boundary points, since 
the boundary is now a curve enclosing a two dimensional space. Let the spatial coordi- 
nates be (x, y), in a general curvilinear coordinate system (not necessarily rectangular). 
The solution U tJ is obtained at the points (x,-,y,) on a rectangular grid with equal spac- 
ings (Ax, Ay) between successive points in each direction. Interior points have 0 < * < I, 
® — 3 — J * The boundaries form a rectangle in the xy plane, and each side of the rectan- 
gle consists of one or more layers of boundary points. The boundary surfaces intersect at 
four corners, each of which consists of one or more corner points. Away from the corners, 
each boundary point has an associated normal and tangential direction (there would be 
two tangential directions in three dimensions), while at the corner points each direction 
is normal. 

The original conservative system of equations is given in (5.26). For definiteness, let 
us consider the y boundaries, defined by the surfaces y=constant, which have the index 
values j < 0 or j > J. Then the x derivative, whi.Qh is in the tangential direction, can be 
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evaluated numerically as an interior term. The y term is in the normal direction, however, 
and must be put in characteristic form so that the appropriate boundary conditions can 
be imposed. Thus we write Eq. (5.26) as 




^T -1 


dU 

MT— + C„ 

ay 




(5.34) 


at the y boundaries. Abbreviating the quantity in parentheses as -3U /dt v , we must 
evaluate d\J/dt v , as given by 


T^ + MT^ + TC f = 0, 
dt v ay 


(5.35) 


to provide boundary conditions for (5.34) at y boundaries. Next define the quantity ,Mfc! 

(5.36) 


M k = 


and compute d\J jdt v from 


for outgoing waves, 


0 


for incoming waves, 


3TJ 

m fc — + M k + m fc C„ = 0. 

ot v 


(5.37) 


The spatial derivative in (5.36) is approximated by the one sided difference formulas 

1 


dV 

dy 


tj ~ A U,j) 3 <0, 

= i>J- 


(5.38) 

(5.39) 


Given d\3 jdt y , we compute dTJ/dt from 


™ ?I + C r = p^ 

dt dx 1 dt u 


(5.40) 


At the x boundaries the y derivatives are evaluated in conservative form by centered 
difference approximations, while the x direction terms are put in characteristic form as 
above. At the corners both directions are normal, and all terms are put in characteristic 
form, as in (5.33). 

5.5 Numerical Solution of the Interior Problem 

The problems considered in this chapter are of two types. The first is the one dimensional 
fluid dynamics problem, in either rectangular or spherical coordinates. The second is a 
two dimensional problem in rectangular coordinates. In both cases the solutions may 
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be discontinuous, and it is necessary to add dissipative terms to the finite difference 
approximations in order to damp nonphysical oscillations around the discontinuities, as 
discussed in Chapters 2 and 4. 

The one dimensional system, in conservative form and with dissipative terms included, 
may be written 


dm 

dt 


dp 

dt 


1 d , _ . 

H —rpt) = 

r n dr K ’ 


Id. „ , 

+ + 

r n dr 


dp 

~&r 


dt 

dt 


+ ^[ rn ( e + P)H 



K 

r* 


m, 


(5.41) 

(5.42) 

(5.43) 


where p is the density, v the velocity, m the momentum density (m = pv ), e the energy 
density (e = + p/(7 — l)), and p the pressure of the fluid. The equation of state can 

be written 


P = 


h -!)(«- \pv 2 ) > 


(5.44) 


where 7 is the constant ratio of specific heats. The spatial coordinate is r, and the 
coordinate system is specified by n (n = 0, 1, or 2 for rectangular, cylindrical, or spherical 
coordinates, respectively). The dissipation coefficient re and finite difference methods are 
given in Chapters 2 and 4. Fourth order finite difference approximations are made to 
the spatial derivatives in the convective and pressure gradient terms. The resulting 
semi-discrete equations are integrated in time according to the method of (2.31). The 
boundary equations (5.25) are integrated along with the interior equations, using the 
same time stepping scheme. 

The two dimensional system, also in conservative form and containing dissipative 
terms, but in rectangular coordinates, is 


dp 


+ x ~{pv x ) + ir-ipv,,) 


dt dx 


dy 


±( V 

dx V dx 


+ 


±( dp\ 

dy \ K dy ) ’ 


dm x d d dp d ( dm x \ d ( dm x \ 

+ 3=K«J + YyM + 51 = a; Hr) + a » ) ’ 


dt dx 


dx 


dm v d d 

dt + dx^ mvVx) + d^ {mv 


, dp d ( dm v \ d ( dm v \ 

”“ ) + ai; = a;( K ^rj + ar(' c Vj' 


(5.45) 

(5.46) 

(5.47) 


de d d . d ( de\ d ( de\ , , 

dt + dx^ + + dy [(C + ~ dx (/"dx] + dy \ K d^ ) ’ (5,48) 

where (m x , m v ) and (v x , v v ) are the momentum density and velocity vectors, respectively. 
The equation of state is 


P=(T- 1) e - y (y\ + vfj . 


(5.49) 
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5.6 Characteristic Equations for Fluid Dynamics 

The boundary conditions require that the fluid equations be put in characteristic form 
at the boundaries, so we begin the boundary specification for fluid dynamics problems 
by finding the characteristic form for the fluid equations. 

In the one dimensional case we can write the fluid equations in the form of Eq. (5.2) 
with 


U 


( \ 

P 


\ ° ; 


( 


where s is a measure of the entropy 


and c is the speed of sound 


v p 0 

d v JL 

P P> 

0 0 v 


VP 7 , 


c = 


~ r pv 

0 

0 


A 


„2 _ 


tp/p- 


(The dissipative terms are set to zero at the boundaries.) 
The eigenvalues of A are 


(5.50) 


(5.51) 

(5.52) 


Aj = t; — c, A 2 = v, A 3 = v + c, 


(5.53) 


and the left eigenvectors are 

1, = (-c, p, -£) , 1, = (0, 0, 1), I, = (c, p, i) . (5.54) 

Taking s as a primitive variable simplifies the eigenvalue calculation, but is incon- 
venient for numerical work. Therefore we eliminate s in favor of p and p and get the 
characteristic equations 


dp dv 
~dt~ PC dt 


+ Ai 



n , 

H — pc 1 v = 0 , 

r 


(5.55) 



where the new primitive variable vector TJ has components p , v, and p. 
Given dXJ/dt at the boundaries, <9U jdt is obtained from (5.3) by 

dp dp 

~dt ~ ~di' 


dm 

~dt 


dp dv 

v Tt +p ^ 


(5.56) 

(5.57) 


(5.58) 


(5.59) 
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dt 1 2 dp dv 1 dp 

— — = — ir ^ pv 1 . 

dt 2 dt dt 7 — 1 dt 


(5.60) 


In two dimensions the fluid equations may be written in the form of (5.33) with 


v x p 0 0 

4 v* 0 


0 0 v~ 0 


0 0 0 v. 


»B = 


Vy 0 fi 0 

0 Vy 0 0 

^ 0 Vy *- 

p • p* 

0 0 0 v„ 


= C„ = 0. (5.61) 


The eigenvalues of A and B are 


Ai — v x c, A2 — A3 — v x , A4 — v x + c, 


(5.62) 


Pi ~ Vy — C, P2 = Hz — Vy, p 4 = Vy + C. 

The y direction characteristic terms, in the form of (5.35), are 


(5.63) 


-pc— + pi 


dp dvy 


(5.64) 


+ P2^ = 0 , 


(5.65) 


2 d P , „ ( dp 2 dp 


dy dy 


(5.66) 


An analogous set holds for the x direction. 

The conservative and primitive time derivatives are related by 

dp dp 
~dt = ~dt' 


(5.67) 


(5.68) 


i x dp dv x 
t “ Vx dt + P dt ’ 


dp dvy 


- o ( v l + + P 


1 dp 
7 — 1 dt ' 


(5.69) 


(5.70) 


(5.71) 
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5.7 Boundary Conditions for Fluid Dynamics 


In the one dimensional case, we solve Eqs. (5.41)— (5.43) in the interior using centered 
finite difference approximations for the spatial derivatives, and an explicit ordinary dif- 
ferential equation solver to integrate the time derivatives of the conservative variables. 
The interior algorithm requires data at the boundary points, which are obtained by solv- 
ing the combined characteristic and nonreflecting equations at those points. The details 
of the boundary calculations are given below. 

We first write the boundary equations as 


dpi 

dt 


dvi „ n * „ 

— PiCi~J7 + -Cl i + ~Pi c i v i ~ 0, 

dt Tf 


(5.72) 


dpi 

dt 


2 dpi 




= 0 , 


dpi avi 

dt + piCi dt 


+ Ci i + —Pic, Vi = 

n 


(5.73) 

(5.74) 


where each Zti is set to zero if At at r< is directed inward (nonreflecting condition), or is 
computed according to the characteristic equations if At at r,- is directed outward. 


Z u 

= ( v < - b*+i 

- Pi - PiCi{vi + 1 - v,)] 

i <0, Vi — Ci 

< 

o, 

(5.75) 


II 

'cT 

1 

n 

1 1 

sr 

i 

Pi— 1 - PiCi{Vi ~ Vi-l)] 

i > I, Vi - c, 

> 

0; 

(5.76) 

Ui 

= «.£; [p.+l - P. - 

- c*(/>i+l - /»<)] * <0 » 

«i < 0, 



(5.77) 


= v, Tr [ p ' _ f ’- 1 ‘ 

- c]{pi - Pi-i)] i> I, 

Vi > 0; 



(5.78) 


= fa + Iw+i 

- Pi + PiCi{v i+ 1 - «,•)] 

* < 0, Vi -1- Ci 

< 

o, 

(5.79) 


= («i +Ci)^{pi - 

Pi— 1 + PiCi{Vi ~ t7,-i)] 

i > 7, V{ + Ci 

> 

0. 

(5.80) 

Given the 

Z values, the time derivatives of the primitive variables are 




II 

"T3 1 ^ 

— ~(Zsi + Zu ) Picf 

Z T% 

Vi, 



(5.81) 


dvi 

dt 

- £«). 

ZpiCi 




(5.82) 


II 

7} + £s ‘) ’ 




(5.83) 


and Eqs. (5.59) and (5.60) provide drrii/dt and dtijdt at the boundaries, to be integrated 
in time along with the interior values. 

The two dimensional case is similar. We solve Eqs. (5.45)-(5.48) in the interior, 
using centered finite difference approximations for the spatial derivative, and an explicit 
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integration method for the time derivatives. We now have four boundaries, defined by 
x = x min > x = imoxi y = 1/mmj and y = y max . Since all four boundaries are treated in a 
similar fashion, it will be sufficient to look at the y=constant boundaries. 

We begin by writing the fluid equations at the y boundary points as 


-x d 


3m. 
dt dx 


dp , d . . dp d ( dp\ 

(rn , „ \ , dP dm x d f 
dm v d dm v d ( dm u \ 

-ar + Tx = 

de d de d ( de\ 

57 +^|(e + pK]- — = 


(5.84) 

(5.85) 

(5.86) 

(5.87) 


where the d\3 / dty terms are the contributions to dX3 jdt due to derivatives in the normal 
(y) direction. The x derivatives are in the tangential direction, relative to the boundary, 
and are evaluated just as in the interior. 

We need to compute the dXJ /dty terms in Eqs. (5.84)— (5.87), and begin by writing 
Eqs. (5.64)-(5.67) in finite difference form, as 


dpn 

dt u 


P*j c ij 1*7 — 

(5.88) 


(5.89) 


(5.90) 

dvuff 

"1“ Pij c tj ^ — 0, 

(5.91) 


dt 

dpij 

dt u „ 

where each Mkij is set to zero if the local y direction characteristic velocity, /ij., is directed 
inward (nonreflecting condition) , or is computed according to the characteristic forms 
below if pic is directed outward: 


•Miu - i v vij ~ [P.y+i - Pij ~ PijCij{vyi j+1 - tv,)] j < 0, v vij - dj < 0, (5.92) 


- (Vyij dj ) ^ [ Pij - Pij -I - Pijdj (vyij - Wyijf-i)] j > J, Vyij - dj > 0 ;( 5 . 93 ) 

M 2i j = Vyij-±-(v xi j +1 - v xi j) j < 0 , v^j < 0 , ( 5 . 94 ) 

1 . 

- Vxij- 1) 3 > J, Vyij > 0; (5.95) 

Msij = Vyij — [p ij+l - Pi j - c?j(pij +l - />,,)] j < o, Vyij < 0, (5.96) 
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= u v ij” [p,y - Pij- 1 - {Pij ~ Pij-i) ] 3 > v v»'i > (5.97) 

Miij = ( v vij + [po+1 - P.J + Pi) c ij{ v vij+ 1 - 3 < °> U V*'J + c »'i < °> ( 5 - 98 ) 

= ( v v*j + c ‘>)^ fotf - Pi}- 1 + ~ v v*i-i)l 3 > J, v vH + c <> > °-(5.Q9) 

Given the M values, we compute dV/dt v from 


dpi 

dt v 

dv 
dt j 

dv 


'•} _ 


V»J 


XIJ 


dt y 

dpij 

dt*! 


Finally, dm x ij/dt v , dm. 


1 

2 PijCij 


* ^ , ¥ , y*j iv' -* , * * * ri f \ fr o*r\ A. 

their values are used in the finite difference approximations to Eqs. (5.84)— (5. 87j . A 
similar process is followed at the x boundaries. 


+ •Ml*;) 5 

(5.100) 

( -M “ «Mli;)> 

(5.101) 


(5.102) 

1 + Mso) • 

(5.103) 

y are calculated 

using (5.69)-(5.71), and 


5.8 Test Problems 

It is instructive to see how well the combined interior and boundary methods described 
up to this point work in practice. To do so, we consider some test problems, in one and 
two dimensions. 

As in Chapter 4, the one dimensional test problems are evaluated on a grid of unit 
length, divided into 100 subintervals ( I = 100, A r = 0.01). The value of 7 in the equation 
of state is 5/3 throughout. The figures show the analytical solutions (solid lines) and 
numerical solutions (dots) for the density p, the pressure p, the momentum density tn, 
and the velocity v, at a time t. The Courant number <x, 

o— h dimension), 

A r x (5.104) 

+ ^3L j (2 dimensions), 

Ax Ay J 

is set to 1 throughout, where c x = max(|v z | + c) over the grid (and similarly for c v and 

*). 

The first problem is a single shock wave moving into a uniform stationary medium. 
The shock starts at x = 0.5 and has a positive velocity. The pre-shock state has p = p = 1 
and v = 0. The problem is determined by the pre- and post-shock values of p, p, and 
v, and the shock velocity V*, and must satisfy the three shock jump conditions [12]. 
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Thus we have one more free parameter, chosen to be the Mach number M = v/c of the 
post-shock flow. The Mach number is related to the ratio R = p,/p u of the shocked to 
unshocked pressures by 


— {R ~ 1} 

7.R7 + 1 + (7 - 1)/?' 


(5.105) 


As R -» oo, M 2 —> M£ oi = 2/[ 7 (7 - 1)], where M max = 1.3416 for 7 = 5/3. 

A reflection is generated when an outgoing subsonic shock wave crosses a boundary, 
as demonstrated by Hedstrom [17]. The reflection takes the form of a constant amplitude 
perturbation to the post-shock solution, which travels inward from the boundary at the 
speed of sound relative to the moving fluid. A convenient measure of the reflection is the 
relative error in the pressure after the shock has crossed the boundary, defined as the 
difference between the numerical and analytical pressure values, divided by the analytical 
pressure. Table 5.1 gives the pressure ratio R and reflections as a function of the Mach 
number, for subsonic shock waves modeled with the dissipation specified by k = 0.35. (A 
significantly smaller value of k results in large oscillations near the shock, while a larger 
value spreads the shock jump out over many grid points.) 

The reflections are small, nowhere exceeding 1%. The worst case is the M = 0.98 
shock, with a reflection of 0.82%. It is interesting to note that the reflection decreases 
as Af increases from 0.99 to 1.0, although the shock jump increases. 

Hedstrom [17] observed a reflection of 12% in the velocity profile of his Figure 6. It is 
not certain why there is such a large discrepancy between his results and those presented 
here, but a likely culprit is the mismatch between his interior and boundary methods. 
He used the Lax-^A^end^off method for the interior. At the boundary, he used first order 
approximations to the spatial derivatives, and 


= Si k ' +1 “ v?) + ° (A ‘> < 5 - 106 ) 

for the time derivatives. The two time integration methods are quite different. In con- 
trast, the four step time integration scheme of (2.31) was used for both interior and 
boundary points here. 

When the flow behind the shock is supersonic (M > 1) all characteristics point to 
the right, and no signals can propagate to the left. Thus no reflections can be produced, 
and none are observed. 

The next example is the shock tube problem, as described in section 4.9. At time 
t = 0 the system consists of two spatially constant, stationary states, adjoining at x = 0. 
The left state (x < 0) has p = p = 1, while the right state has p = 0.125, p = 0.1. 

Figure 5.1 shows the numerical solution at step 100 and t = 0.438, with k = 0.3. The 
rarefaction wave has passed through the left boundary, and the shock wave has passed 
through the right boundary. No reflections are visible on either side. The boundaries are 
well behaved. 

The next problem considered is the Sedov solution for a spherically symmetric explo- 
sion, also described in section 4.9. Figure 5.2 shows the explosion at step 1000, t = 3.069, 
with k = 0.35, and with the vertical scales magnified relative to the example in Figure 4.8. 
The shock wave has passed out of the domain, and the velocity at the boundary has de- 








Figure 5.1: Shock tube at step 100 
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creased from supersonic (immediately behind the shock) to subsonic. A perturbation 
has developed at the boundary and is propagating inward, as can be seen most clearly 
from the velocity curve. The discrepancy between the numerical and analytic solutions 
is presumably due to the fact that the similarity solution for the explosion really does 
contain an inward propagating wave, which is suppressed by the boundary conditions. 
The resulting numerical solution is a valid one, but it is not the solution to the explosion 
problem at late times. 

The impossibility of properly specifying boundary conditions for all problems with the 
nonreflecting prescription is further illustrated by the following problem, the homologous 
expansion of a uniform medium. At time t = 0 the density and pressure are uniform, with 
p = p=l. The velocity is linear, with u = x/t 0 , to chosen to be 1. The region studied is 
-0.5 < x < 0.5. The density and pressure decrease with time but remain uniform, while 
the velocity also decreases but remains linear in x. The flow at the boundaries is always 
subsonic and directed outward for this set of initial conditions. 

The problem has a simple analytical solution, given by 


°( 1 + h) ' 

(5.107) 

( t\-'> 


} ( 1+ ^) ■ 

(5.108) 

X 


t + io 

(5.109) 


One can verify by direct substitution that this solution does not satisfy the nonreflecting 
boundary conditions. For example, at the right boundary (x = b) the flow is subsonic 
and directed outward. The two characteristic equations representing outgoing waves 
(Eqs. (5.56) and (5.57)) hold as written, while the absence of an incoming wave is imposed 
by the nonreflecting condition of Eq. (5.72) with = 0, which can be written 


dp 

dt 


dv 

dt 


pc— = 0. 


The analytic solution does not satisfy the nonreflecting condition, so the nonreflecting 
condition will not produce the desired numerical solution. 

Figure 5.3 shows the expansion problem at step 50, t = 0.305, with jfc = 0. The 
numerical and analytical solutions diverge markedly near the boundaries, although they 
match well in the middle portion. The discrepancy grows with time until the two solutions 
disagree everywhere. 

We can use information about this particular problem to specify better boundary 
conditions. In particular, the pressure gradient is zero everywhere, and the velocity 
satisfies 

dv dv 

di + v ai = 0 ’ ( 51 io) 

which is in characteristic form and describes outflow at the boundary. If (5.110) is used 
in place of the nonreflecting boundary condition, we have three characteristic equations 


rionENTun 
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describing outgoing waves. (By throwing out the pressure gradient term, we have ex- 
cluded sound waves from the problem. The only characteristic velocity left is the fluid 
velocity. The evolution of the velocity profile then determines the rest of the solution). 
The resulting numerical solution matches the analytical solution everywhere. 

This result is not meant as an endorsement of Eq. (5.110) as a boundary condition 
in general (it does nothing useful for the explosion problem), but simply illustrates that 
nonreflecting boundary conditions cannot be expected to give the desired results on 
problems which do contain incoming waves. Information specific to such problems may 
be used to produce more useful boundary conditions. 

The final test problem is two dimensional. It is a planar Newtonian shock, traveling 
in rectangular geometry. The grid has Ax = Ay = 1/130. The shock is traveling toward 
the upper right, at a 45 angle with respect to the x axis. The initial distance Rq between 
the shock front and the origin is 0.8. The unshocked density and pressure are p = 1, 
P with v x = v v — 0. The Mach number of the flow behind the shock is M = 0.95, 
picked to roughly maximize reflection errors. The dissipation used is k = 0.3. 

Contour plots of the density, pressure, and velocity fields in Figures 5.4, 5.5, and 5.6 
show the time evolution of the solution. The arrows in the velocity plot show the direction 
of the flow, and have a length proportional to the magnitude of the velocity. The ratio 
of successive contour values is 1.01. The figures are at step 0, t = 0; step 100, t = 0.098; 
and step 200, t = 0.197, respectively. Step 0 shows the initial conditions. Step 100 
shows the shock shortly before it reaches the corner. Small boundary perturbations can 
be seen near the edges of the shock. Step 200 shows the solution shortly after the shock 
has left the grid. The reflection from the corner has propagated inward, and at its peak 
amounts to about 6% of the post-shock pressure profile. The reflection is greater than in 
the one dimensional case, perhaps because the corner is subject to reflections from two 
coordinate directions, but not enough to obscure significant features of the post-shock 
flow (if there were any). 


Y FIX I S 


5.8. TEST PROBLEMS 


91 



Figure 5.4: Angled shock at step 0 
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Figure 5.5: Angled shock 
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Figure 5.6: Angled shock at step 200 
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5.9 Exercises 

1. Derive the primitive form of the fluid equations given in Eqs. (5.50) from the 
conservative form, in one dimension and in rectangular coordinates. 

2. Starting with Eqs. (5.50), assume isentropic flow (s = constant), and show that 
only two characteristic equations result. What are the characteristic velocities? 

3. Express the results of problem 2 in the form of Eq. (5.15). The two wave amplitudes 
Vi and V 2 are called Riemann invariants in fluid dynamics. Each Riemann invariant 
is constant along its corresponding trajectory (v - c or v + c) in the xt plane. 
Riemann invariants do not exist if the flow is non-isentropic. 

4. Find the right eigenvectors of matrix A in Eq. (5.49). Form the similarity transfor- 
mation matrices S and S -1 , and verify the similarity transformation of Eq. (5.11). 

5. Derive the characteristic velocity components for the two dimensional fluid problem 
of Eq. (5.61). 

Perform the shock wave calculation described in section 5.8, for a Mach number 
M = 0.9. Do you see a reflection? If so, how big is it? 

7. Verify that the homologous expansion solution of Eq. (5.107)-(5.109) does not 
satisfy the nonreflecting boundary condition at the right boundary. 


Appendix A 

The Equations of Fluid Dynamics 


Chapter 4 described the basic concepts of tensor calculus, the tensor formulation of the 
fluid dynamics equations, and artificial viscosity. This appendix provides a complete 
list of the relevant equations for fluid dynamics in rectangular, cylindrical, and spherical 
coordinate systems. 

A.l Coordinate Geometries and ^Metric Tensors 

The geometry of the coordinate system is described by the metric tensor g a b. The metric 
tensor is a symmetric second rank tensor, whose components are usually functions of 
the coordinates. The metric tensor defines the physical distance ds between two closely 
space points, at x° = {x\x 2 ,x 3 ) and x“ + dx a = (x l + dx\x 2 + dx 2 ,x 3 + dx 3 ), in terms 
of the coordinate differentials dx a . This distance is given by given by 

ds 2 = g ah dx a dx\ (A.l) 

The symmetry of the metric tensor implies that, in n dimensions, only n(n + l)/2 
of its n 2 components may be unique (six components in three dimensions). A general 
metric in three dimensions may be written as a 3 x 3 array of numbers, 

<7ll 9 12 ffl3 
921 922 923 

931 932 933 

where by symmetry 

gab = 9ba • (A- 3 ) 

In orthogonal coordinate systems, where the coordinate axes are perpendicular to 

C~ & 
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each other, the metric tensor is diagonal, and may be written 


9ab = 


h \ 


h\ 


hi 


(A.4) 


where h u h 2 , and h 3 are the scale factors for the coordinate directions (and the superscript 
2 represents the square of the number, not a tensor index). 

The covariant derivative of a tensor in coordinate direction a is denoted by the sub- 
script ; a. The following are useful examples: 

** 

£ H 

II 

o 

(A.5) 

B a . — ■ P“ DC 

il “ dx> + r “ B • 

(A.6) 

B _ dBa p. „ 

tia ' b dx b 1 ai ^ e ’ 

(A.7) 

dC ab 

c<lh , e — dxe +r3 e c* + r5,c* 

(A.8) 

c ‘ >x = dJ + r * cd ‘ ~ r ^ c “<'> 

(A.9) 

c.,, = ^ - Tifi* - Tic. d . 

(A. 10) 

The connection coefficients T bcy used to evaluate covariant derivatives, 

are given by 


r? c = g^rac, r 


dbe 


■K 


dg<tb dgdc _ dgbc 

dx c dx b dx d 


)■ 


(A.11) 


and are symmetric on the last two indices (FJ. = r® 6 ), due to the symmetry of the metric. 
Note that T bc is not itself a tensor. The connection coefficients also satisfy the equation 


ta gdx b dx b ln9 ' 


(A.12) 


where g = 0det^jf. 
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A. 2 The Coordinate Invariant Fluid Equations 

The density (p), pressure (p), and total energy density (e) of a fluid are scalar fields and 
are invariant under coordinate transformations. The coordinate velocity (u“), momentum 
density (m a = pu a = pg a bU b ), mass flux (pu°), and total energy flux ([e + p]u a ) are vector 
(first rank tensor) fields. The momentum flux (pu a u b + ptf) is a second rank tensor. 

The coordinate velocity u° = dx a Jdt is the rate of change of a fluid element coordinate 
z° with time, and is generally not the same as the physical velocity v a , which is the rate of 
change of distance along coordinate axis z°. In an orthogonal coordinate system, where 
the coordinate axes are at right angles to each other, the metric g a b is diagonal, with 
diagonal elements g aa = h\. In such a system the physical and coordinate velocities are 

related by v a = h a u a (not summed). _ # 

Closing the system of equations requires auxiliary algebraic equations as well as the 
differential equations given below. The relationship between momentum density and 
velocity is one such equation: 

u a = g ab Ub = g ab m b /p. (A. 13) 

Also useful is the velocity magnitude v, given by 

u 2 = u„u° = g ab u a u b = v 2 + Uj + V 3 . (A. 14) 

The equation of state connects the pressure p to the thermal energy density e according 

to 1 1 

p = (7 - l)c, e = e - - pv 2 = e - -m a u°. (A.15) 

In simple perfect gas problems 7 is the ratio of specific heats, and is a constant (equal 
to 5/3 for a monatomic gas). In more complicated problems (such as relativistic fluid 
dynamics) 7 may be a function of the local fluid state, and is not the ratio of specific 

heats 

The equations given in this appendix include dissipative terms for use as artificial 
viscosity in finite difference calculations. The dissipation coefficient, k, is taken to be 
nonzero only because the grid resolution used in the calculations is finite. To get the 
dissipation-free equations, simply set k = 0 throughout. The definition for k is given m 
Chapter 4 and will not be repeated here. 

The coordinate invariant fluid equations, with dissipation, are 


+ 0™ a );« = ( K P' a );a’ 

(A.16) 

— + [(e + p)u a ];a = (kO ; o , 

(A. 17) 

^ + (m a u 6 + p6 b a )-b = (*rri b ) b , 

(A. 18) 

in terms of covariant derivatives, or 


dp Id, n , 1 d ( ab d P\ 

dl + g Ss« (9 ' m ) ~ 9 dx‘ l 9 * 5 d* 1 )' 

(A.19) 
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l + ^ w ' +p)u * 1 = ]£■ (‘"'*“ 1 ?) - 


(A. 20) 


dm a 

dt 


1 d b dp 1 dg ei 

+ 9 m » U + ~sT~a + o P U C U d~^— 

g dx b dx a 2 dx° 


1 3 


g dx c 


(I? - r -‘^)] - r ~ [»"* (I? - >*»*) 


, (A. 21) 


where g = y'l det(y at )|. 

The following sections give the fluid equations explicitly in rectangular, cylindrical, 
and spherical coordinates. It should be noted that the forms of the equations given are 
not exactly those obtained from the above calculations. Care has been taken to remove 
terms which cancel each other (especially when the terms individually become infinite, 
at singularities) and to cast the equations in a form well suited for finite difference 
calculations. The direct translation of the tensor forms given above into finite difference 
forms leads to serious cancellation errors. One should always eliminate cancellations 
and singularities by hand to the greatest degree possible before resorting to numerical 
approximations. - 


A. 3 Rectangular Coordinates 

The coordinates are x a = (x,y,z). The metric and its inverse are 





\ 


l 


1 

Qab — 

l 

ab 

y 9 = 

1 


i * J 


l 1 1 


(A. 22) 


and the metric factor g — 1 . The contravariant (coordinate), covariant, and physical 
velocities are the same in this (and only this) coordinate system: 


u° = (i,y,z), 

u a = (x,y,z), (A. 23) 

v a = {i,y,z) = (u*,U y ,U*). 


A. 3.1 Connection Coefficients 

All connection coefficients are zero: 

K = o. 


(A.24) 


AA. CYLINDRICAL COORDINATES 
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A. 3. 2 The Fluid Equations with Dissipation 

Set k = 0 throughout for the inviscid equations. 


Density 




k ' ' V ' ' dz ^ 7 

d ( dp\ ^ d ( dp\ ^d_( dp\ 

dx V dx) dy \ dy) dz \ dz) ’ 


Energy 


ft + £l( e +ri“’1 + |;[( e+ p> u,, l A l(e + p)u ‘ 

d ( de\ d ( de\ d ( de\ 

= te{ K te) + ai{ K ai) + a-z{ K T*)’ 


x momentum 


dm x , & i x \ / v\ , & , z\ , dP 

■aT + ei (m ’ u > + ^< m * u > + aJ (m *“ 5 + ai 


d ( dm x \ d ( < 


dy ) dz \ dz 


y momentum 


dm v d . d . V \ , d . dp 

~aT + > + ^ (m ' K 1 + ai {m ‘ u J + ai 

d ( dm v \ d ( dm v \ d ( dm y \ 

Tx (^"aT ) ~dy v K_ ^T J ~dz J 


z momentum 


dm z d d . d . dp 

~aT + sJ (m - u > + ^ (m *“ > + ai (m -“ > + 

d ( dm z \ d ( dm z \ d f dm z \ 

~ dx \ dx ) + dy \ dy ) dz \ dz ) ’ 

A. 4 Cylindrical Coordinates 

The coordinates are x a — (r, <j>,z), and are related to rectangular coordinates by 


x = rcos<j>, r = yx 2 + y 2 , 
y = r sin <f>, tan <i> = y/x, 
z = z, z = z 
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The metric and its inverse are 



A 


\ 


1 


1 

Qab — 

r 2 

_afc 

» 9 = 

r- 2 






(A. 31) 


and the metric factor g = r. The contravariant (coordinate), covariant, and physical 
velocities are related by 

« a = (^,<M), 

u a = (r,r 2 <£,z), (A. 32) 

v a = (f,r<£,z) = (u r ,ru*,u z ). 


A.4.1 Connection Coefficients 

The connection coefficients are 



All other r£ c = 0. 


(A. 33) 


A. 4. 2 The Fluid Equations with Dissipation 

Set k — 0 throughout for the inviscid equations. 


Density 


dp 

dt 


+ 


Id. .. d . d . .. 

7a7 (r/,u) + 5J ( ' >u ) + ai ( ' ,u ) 
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Energy 
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\ i__a / am,\ a / am,\ (A . 38) 

) r*d<f>\ d<f> J dz\ dz J ' 


A. 5 Spherical Coordinates 


The coordinates are x a = ( r,6,<j> ), and are related to rectangular coordinates by 


x = r sin0 cos<£, 

y = rs\n6s\n<f>, 
z = r cos 6, 


r = \jx 2 + y 2 + z 2 , 

tanS = V£!±Z 

z 

tan <}> — yjx 


(A. 39) 


The metric and its inverse are 



r i 

\ 


f i 

\ 

9ab — 

r 2 


ab 

, g = 

r~ 2 




r 2 sin 2 6 J 



r -2 sin -2 0 , 


(A. 40) 


and the metric factor y = r 2 sin 0 . The contravariant (coordinate), covariant, and physical 
velocities are related by 


u a = (r,0, <t>), 

u a = (r, r 2 0, r 2 sin 2 6<f>), 

v a — {r,rO,r sin 0<f>) = (u r ,ru\rsmOu*). 


(A.41) 
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APPENDIX A. THE EQUATIONS OF FLUID DYNAMICS 


A. 5.1 Connection Coefficients 

The connection coefficients are 


r 

r 
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00 = ~ r > 


= -rsin 2 d. 
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r$ ~ 1 Or “ 
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r 

= — sin d cos d 

ii 

•n 

II 

i 

r 

= cot d. 


All other rj c = 0. 

A. 5. 2 The Fluid Equations with Dissipation 

Set k = 0 throughout for the inviscid equations. 
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, (A. 43) 


(A. 44) 


(A. 45) 
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